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I.  DIRECTOR’S  OVERVIEW 


This  report  represents  the  eighteenth  annual  summary  of  The  Ohio  State  University  Joint 
Services  Electronics  Program  (JSEP).  Prof.  Kouyoumjian  has  been  elected  to  the  National 
Academy  of  Engineering.  His  citation  denoted  the  development  of  the  Uniform  Theory  of 
Diffraction.  This  research  was  supported  by  JSEP. 

There  have  been  a  total  of  39  Ph.D.  and  29  M.Sc.  degrees  in  Electrical  Engineering 
obtained  under  partial  JSEP  sponsorship.  There  are  currently  8  Ph.D.  and  3  M.Sc.  students 
being  partially  supported  under  JSEP. 

As  may  be  seen  in  the  Annual  Report  Appendix,  14  reprints  have  been  included  for 
the  period  September  1994  to  September  1995.  In  addition,  10  papers  have  already  been 
accepted  for  publication  in  the  coming  year,  an  additional  6  papers  have  been  submitted, 
and  finally,  16  papers  are  in  preparation. 

II.  DESCRIPTION  OF  SPECIAL  ACCOMPLISH¬ 
MENTS  AND  TECHNOLOGY  TRANSITION 

The  transfer  of  the  compact  range  and  target  identification  technology  initiated  under  JSEP 
support  for  time  domain  studies  continues  to  make  large  advances.  The  large  compact  reflec¬ 
tor  and  a  modern  radar  of  our  design  at  ASD,  Wright  Patterson  Air  Force  Base  (WPAFB) 
has  been  the  subject  of  our  updated  evaluation  techniques  resulting  in  substantial  improve¬ 
ment  in  its  operation,  especially  for  lower  frequencies.  WPAFB  is  interested  in  operating  at 
much  lower  frequencies  than  originally  intended  which  has  required  a  complete  re-evaluation 
of  their  system  at  lower  frequencies.  This  led  to  significant  stray  signal  errors  that  cannot 
be  solved  by  conventional  methods,  such  as  using  12-foot  absorber.  As  a  result  of  our  the¬ 
oretical  JSEP  studies,  we  have  shown  that  much  thinner  material  can  be  built  and  used  to 
solve  the  WPAFB  problems.  An  even  more  advanced  compact  range  has  been  designed  for 
NASA,  Langley  and  is  expected  to  be  built  this  coming  year.  It  was  specifically  designed  for 
measurements  as  low  as  300  MHz.  This  new  range  will  fit  in  a  room  40'  x  40'  x  80'  instead 
of  the  previously  required  120'  x  120'  x  360'.  The  cost  savings  in  this  case  are  enormous. 


1 


New  absorbers  and  feeds  have  been  designed  which  are  to  be  incorporated  in  the  NASA  and 
WPAFB  ranges. 

We  continue  to  assist  Rockwell  (Tulsa)  to  update  their  RCS  facilities.  Our  latest 
pulsed/CW  radar  system  has  been  installed  at  Rockwell  during  the  previous  year.  This 
work  is  on  a  subcontract  to  the  ESL  from  the  Air  Force.  These  and  other  advances  were 
only  possible  because  of  the  initial  JSEP  support.  This  continues  to  be  a  case  where  a  small 
investment  of  basic  research  funds  have  been  leveraged  to  generate  much  larger  support  and 
have  achieved  major  contributions  for  DoD.  This  has  also  led  to  OSU-ESL  involvement  in 
the  study  of  Ultra  Wide  Band  radar  systems. 

The  Compact  Range  Technology  has  also  resulted  in  an  antenna  for  the  detection  of 
anti-personnel  mines.  This  is  a  small  reflector  antenna  which  has  a  secondary  focus  in  the 
vicinity  of  the  ground.  The  extent  of  the  focussed  spot  has  been  evaluated  using  the  Gaussian 
Beam  representation  (see  Diffraction  Work  Unit)  for  the  fields  of  a  reflector  antenna  and  the 
analysis  shows  the  spot  size  to  be  maintained  for  a  substantial  region.  This  is  clearly  the 
most  successful  antenna  found  for  this  task.  The  results  obtained  in  the  Finite  Element 
Work  Unit  have  also  been  used  to  evaluate  the  scattered  fields  from  typical  anti-personnel 
mines.  We  are  now  seeking  techniques  to  be  used  to  identify  these  mines. 

One  of  the  prime  target  identification  candidates  is  the  Complex  Natural  Resonances 
(CNR’s)  which  was  also  funded  by  JSEP  a  decade  or  so  ago.  Such  CNR’s  have  proven  useful 
in  identifying  UneXploded  Ordnance  (UXO’s)  in  the  past  year.  A  demonstration  was  held 
at  Tyndal  Air  Force  Base  over  a  series  of  UXO’s  buried  in  sand.  They  were  easily  detected 
by  the  Ground  Penetrating  Radar  and  their  lengths  were  correctly  deduced  by  the  CNR 
identification  techniques. 

The  UTD  solutions  which  have  been  developed  for  analyzing  the  phenomena  of  EM 
radiation/scattering  within  paraxial  regions  of  smooth  elongated  convex  bodies,  and  also 
those  UTD  solutions  which  are  currently  under  development  for  extending  these  paraxial 
solutions  to  include  end  caps  in  otherwise  smooth  elongated  convex  bodies,  are  of  significant 
importance  in  predicting  the  effects  of  the  shadowing  of  a  jet  engine  inlet  by  the  fuselage  of 
an  aircraft,  when  it  is  illuminated  near  nose  on  by  a  radar.  It  appears  from  our  preliminary 
analysis  that  even  though  a  significant  part  of  an  inlet  is  in  the  shadow  of  the  fuselage  to 
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which  it  is  attached,  the  indirect  illumination  of  the  shadowed  part  of  the  inlet  opening  that 
can  result  from  the  diffraction  within  paraxial  region  of  the  fuselage  is  almost  as  strong  as 
the  field  that  would  have  directly  illuminated  the  inlet  cavity  opening  if  the  fuselage  was 
absent.  This  work  is  presently  being  conducted  for  the  Air  Force. 

Work  is  still  continuing  on  the  development  of  naval  shipboard  antenna  models  based 
on  the  generalized  ray  expansion  (GRE)  concept  for  the  Navy.  Also  OSU  is  involved  with 
assisting  the  NAVY  by  incorporating  UTD  concepts  into  a  ray  based  code  in  order  to  improve 
its  capabilities  and  efficiency  for  predicting  the  effects  the  topside  environment  on  shipboard 
antenna  performance. 

The  UTD  and  hybrid  developments  leading  to  systematic  approaches  for  controlling  un¬ 
desired  edge  effects  in  the  radiation  by  antennas,  via  material  treatment,  is  continuing  to  be 
utilized.  In  particular,  different  types  of  new  materials  are  being  systematically  incorporated 
into  the  design  for  controlling  the  seriously  detrimental  edge  diffraction  effects  of  platforms 
for  microstrip  antenna  elements  in  GPS  applications.  This  work  is  being  performed  for 
Trimble.  Resistive  cards  are  also  being  incorporated  in  a  novel  fashion  into  the  design  of 
reflector  antennas  for  controlling  edge  diffraction  effects  without  modifying  the  feed  illumi¬ 
nation.  This  work  is  being  done  for  Rockwell.  In  addition,  R-cards  are  being  employed  for 
making  measurements  on  test  antennas  mounted  in  a  small  size  ground  plane  such  that  the 
edge  effects  of  these  finite  ground  planes  can  be  essentially  removed,  simulating  the  behavior 
of  such  antennas  on  large  bodies  such  as  an  aircraft  fuselage. 

The  method  of  moments  (MM)  solution  for  a  wire  attached  to  an  arbitrary  body  is  be¬ 
ing  integrated  into  the  “Helicopter  Antenna  Radiation  Prediction”  (HARP)  code,  currently 
being  developed  under  NASA  support.  HARP,  which  features  a  sophisticated  graphical  in¬ 
terface,  is  a  collection  of  high  and  low  frequency  computer  codes  for  analyzing  the  radiation 
from  antennas  on  complex  bodies  such  as  helicopters,  aircraft,  autmobiles,  etc.  HARP’s  low 
frequency  code  is  the  “Electromagnetic  Surface  Patch  Code”  (ESP),  which  at  present  cannot 
treat  a  wire  attached  to  an  arbitrary  surface.  Thus,  when  HARP  or  ESP  analyzes  a  monopole 
on  a  helicopter,  the  monopole  must  be  replced  by  a  dipole  “floating”  a  short  distance  off 
the  body.  While  this  is  a  reasonable  approximation  for  far  zone  patterns,  it  is  impossible 
to  compute  such  near  zone  quantities  as  input  impedance  and  mutual  coupling.  However, 
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when  the  solution  for  a  wire  attached  to  an  arbitrary  body  is  integrated  into  ESP,  HARP 
will  be  able  to  better  model  antennas  on  complicated  support  structure,  and  in  particular 
compute  input  impedance  and  mutual  coupling. 

These  programs  are  also  in  use  in  a  recently  funded  Technology  Transfer  Program  by 
the  Ford  Motor  Company.  This  program  is  concerned  with  the  design  of  hidden  automobile 
antennas. 

The  Perfectly  Matched  Absorbing  Layer  (PML)  boundary  condition  is  likely  to  become 
the  choice  for  finite  element  and  finite  difference  codes  in  the  near  future  because  of  its  ease 
of  implementation  and  its  good  performance.  For  frequency  domain  finite  element  methods, 
the  anisotropic  PML  is  probably  the  most  popular  choice,  because  many  of  the  currently 
available  finite  element  codes  already  can  handle  anisotropic  materials.  There  are  researchers 
who  have  already  begun  adapting  the  anisotropic  PML  to  their  codes  including  Professor 
John  Volakis  at  Michigan  and  Professor  Steve  Gedney  at  Kentucky.  At  this  point,  the  use  of 
the  anisotropic  PML  in  the  time  domain  is  uncertain;  however  if  it  can  be  developed,  then 
it  offers  possible  advantages  for  time  domain  finite  element  methods. 

The  currently  available  elements  for  frequency  domain  finite  element  methods  are  not 
very  efficient  for  electromagnetic  modeling.  The  work  on  the  solenoidal  elements  under 
JSEP  funding  has  shown  that  these  elements  have  great  potential  for  significantly  improving 
the  performance  and  efficiency  of  finite  element  codes.  The  positive  initial  work  will  likely 
lead  to  funding  from  Motorola  for  developing  finite  element  codes  for  them.  The  movement 
toward  the  development  of  numerical  methods  specifically  for  electromagnetic  modeling  will 
lead  to  a  new  generation  of  high  performance  and  high  efficiency  codes  in  the  near  future. 

It  becomes  clear  from  the  above  that  a  major  portion  of  technology  transition  of  JSEP 
research  takes  the  form  of  additional  supported  research  from  other  DoD  agencies.  Other 
forms  of  transition  include  former  graduate  students,  who  upon  graduation  are  employed  in 
DoD  related  positions,  and  publication  of  JSEP  research  results.  Yet  another  form  of  such 
transition  is  represented  by  computer  codes  that  incorporate  the  results  of  our  research. 
These  complex  codes  are  available  to  DoD  related  industries  for  a  nominal  fee  of  8300.  Last 
year  125  such  codes  were  issued  which  approximately  doubles  the  number  issued  the  year 
before.  This  indicates  a  resurgence  of  interest  in  Electromagnetics  in  the  DoD  community. 
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Our  JSEP  research  continues  to  focus  on  electromagnetics  topics.  The  four  work  units 
included  in  last  year’s  Summary  report  are  each  showing  significant  results. 

The  Diffraction  Work  Unit  includes  several  subtopics.  The  Time  Domain  Uniform  Theory 
of  Diffraction  (TD-UTD)  concept  introduced  last  year  is  showing  excellent  progress.  Recall 
that  the  standard  FD-UTD  was  developed  under  JSEP  and  has  been  very  frequently  used  by 
non-OSU  researchers.  As  noted  in  the  Director’s  Overview,  this  led  to  Prof.  Kouyoumjian’s 
election  to  the  National  Academy  of  Engineering.  This  year,  the  TD-UTD  has  been  extended 
to  include  a  much  more  general  wedge  geometry,  the  slope  diffraction  case  and  the  shadow 
side  solution  for  the  creeping  wave  generated  by  smooth  convex  conducting  shapes.  These 
are  the  very  terms  most  used  by  the  FD-UTD  and  for  a  long  time  formed  the  basis  for  such 
FD-UTD  solutions.  It  is  anticipated  that  the  TD-UTD  will  prove  as  rewarding  for  transient 
solutions  as  the  FD-UTD  has  been  for  steady  state  harmonic  solutions. 

The  concept  of  an  asymptotic  analysis  using  Gaussian  Beams  (GB’s)  when  the  use  of 
rays  becomes  computationally  intensive  or  inaccurate  (i.e.,  at  caustics)  is  continuing.  The 
limitation  on  special  edge  orientation  has  been  replaced  by  one  where  a  curved  edge  can 
be  arbitrarily  located  provided  it  is  locally  plane.  Research  in  this  area  is  continuing  and 
is  expected  to  be  completed  soon.  This  activity  should  lead  to  a  far  more  efficient  and 
physically  appealing  solution  when  applied  to  synthesis  of  antenna  feeds  and  also  synthesis 
of  large  shaped  reflectors  for  generating  multiple,  spot  or  contour  beams.  The  use  of  GB’s 
for  the  treatment  and  design  of  large  planar  arrays  is  to  be  investigated.  They  are  expected 
to  play  a  very  significant  role  in  the  development  of  antenna  synthesis  technique  for  wireless 
communications  systems,  and  the  analysis  of  radomes  used  in  conjunction  with  electrically 
large  antennas. 

The  Integral  Equations  Work  Unit  has  shifted  its  focus  from  the  treatment  of  Artificial 
Media  to  a  novel  time  swept  approach  for  scattering  from  large  bodies.  The  concept  is  to 
sweep  over  the  body  in  finite  segments  with  each  segment  terminated  in  a  tapered  R-card 
to  eliminate  edge  diffraction  terms.  This,  in  essence,  breaks  the  large  body  problem  into  a 
series  of  small  bodies  with  a  substantial  reduction  of  computer  time  and  resources. 

The  Finite  Element  Work  Unit  focussed  attention  on  two  topics  this  past  year  which 
should  have  a  substantial  impact  on  these  powerful  solutions  for  electromagnetic  problems. 
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The  first  is  the  Perfectly  Matched  Absorbing  Layer  (PML)  which  is  used  to  truncate  the 
computation  domain  in  differential  based  numerical  solutions.  The  adaptation  of  anisotropic 
PML  is  considered  for  Time  Domain  Finite  Element  Solution.  This  involves  a  question 
of  dynamic  stability  which  is  examined  carefully.  The  second  topic  being  pursued  is  the 
solenoidal  edge  element.  It  has  been  specifically  designed  to  minimize  errors  due  to  numerical 
dispersion. 

The  Hybrid  Studies  Work  Unit  continues  to  demonstrate  that  it  is  achieving  the  same 
prominence  that  the  original  UTD  has  achieved  by  being  able  to  treat  a  number  of  config¬ 
urations  that  would  otherwise  be  untreatable  by  any  single  form  of  analysis.  This  Hybrid 
Work  Unit,  of  course,  demonstrates  the  synergistic  nature  of  our  JSEP  program. 

As  an  example,  the  design  analyses  and  measurement  of  the  properties  of  an  oscillator 
incorporated  into  an  antenna  has  come  to  fruition  in  the  past  year.  This  incorporates  the 
solid  state  electron  device  and  antenna  technology  as  a  hybrid  solution.  Other  examples 
include  the  use  of  various  forms  of  Green’s  functions  (impedance  loaded  wedges  and  UTD 
solution  of  coated  convex  surfaces)  with  obstacles  for  both  radiation,  radome  and  scattering 
problems.  Another  form  of  Hybrid  Solution  involves  the  Fast  Multiple  Method  (FMM) 
which  is  gaining  increasing  interest  by  the  computational  community.  Modification  of  the 
FMM  makes  use  of  asymptotic  techniques  to  make  this  high  speed  solution  even  faster. 
Several  interesting  computational  techniques,  are  also  being  examined  that  may  ultimately 
be  included  in  the  Hybrid  Library.  These  are:  a)  the  Finite  Volume  Time  Domain  (FVTD) 
Solver  which  could  prove  to  be  superior  to  the  FDTD  Solver  in  some  applications,  and  b) 
a  more  efficient  UTD  transition  function  involving  Fock  type  integrals  that  appear  in  the 
UTD  special  Green’s  functions  which  will  dramatically  speed  up  computations  for  radiators 
mounted  on  convex  material  coated  conducting  surfaces. 
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1.  Introduction 

The  frequency  domain  (FD)  based  uniform  geometrical  theory  of  diffraction  (UTD),  devel¬ 
oped  mostly  at  The  Ohio  State  University  and  primarily  under  JSEP  support,  has  proven 
to  be  an  efficient  analytical  method  for  the  prediction  of  EM  radiation/scattering  from  elec¬ 
trically  large,  relatively  complex  objects  for  which  exact  analytical  solutions  are  virtually 
impossible  to  obtain.  In  addition,  the  UTD  concept  is  simple  and  it  provides  physical  in¬ 
sight  into  the  radiation/scattering  mechanisms.  In  contrast,  brute  force  numerical  solution 
methods  lack  insight  into  the  relevant  physics  and  become  rapidly  inefficient  or  intractable 
as  the  radiating  object  becomes  large  in  terms  of  the  electrical  wavelength. 

Of  course,  the  range  of  phenomena  to  which  the  UTD  may  be  applied  as  an  analy¬ 
sis/design  tool  depends  on  the  availability  of  the  pertinent  canonical  UTD  solutions  from 
which  the  UTD  model  of  a  complex  radiating  object  can  be  “built  up.”  While  several  canon¬ 
ical  UTD  solutions  have  been  obtained  for  perfectly-conducting  boundaries  (e.g.,  curved 
wedges;  smooth  convex  surface;  etc.)  and  recently  also  for  non-conducting  or  material 
boundaries  (e.g.,  edges  in  uniform  material  slabs,  and  wedges  with  impedance  boundary 
conditions,  etc.),  there  still  remain  a  host  of  important  canonical  UTD  configurations  which 
need  to  be  analyzed.  Furthermore,  a  similar  basic  problem  to  that  in  the  frequency  domain 
(FD)  exists  also  in  the  time  domain  (TD),  namely  brute  force  numerical  TD  solutions  cannot 
provide  efficient  or  even  tractable  procedures  for  predicting  the  transient  response  of  short 
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pulse  excited  complex  radiating  objects.  The  need  for  efficient  TD  analysis  based  prediction 
tools  is  motivated  by  the  development  of  short  pulse  or  ultra  wide  band  (UWB)  radars  for 
target  identification,  and  remote  sensing.  Also,  it  is  important  to  predict  the  effects  of  nat¬ 
ural  and  manmade  electromagnetic  pulses  (EMP)  on  EM  radiating  systems.  The  FD  based 
approaches  are  also  cumbersome  for  predicting  transient  response  because  the  FD  solution 
has  to  be  found  at  a  very  large  number  of  discrete  frequencies  and  then  inverted  numerically 
into  the  TD  via  the  IFFT  algorithm. 

Thus,  during  the  last  year,  the  diffraction  studies  unit  has  focused  on  the  two  major  issues 
discussed  above.  In  particular,  we  have  made  significant  progress  in  the  development  of  a  new 
TD-UTD  for  analyzing  early  to  intermediate  time  transient  response  of  some  highly  useful 
canonical  configurations,  and  their  generalizations,  so  that  realistic  and  hence  complex,  pulse 
excited  radiating  objects  can  be  handled  efficiently  and  in  a  physically  appealing  manner 
via  TD  rays.  The  TD-UTD  development  is  in  addition  to  the  important  progress  made  in 
the  development  of  new  and  useful  FD  based  asymptotic  HF  or  UTD  solutions  for  non¬ 
conducting  as  well  as  conducting  objects.  The  increased  use  of  composite  materials  in 
aerospace,  satellite  and  other  antenna  systems  motivates  the  FD-UTD  research  on  material 
effects  on  radiating  objects.  Specifically,  in  the  area  of  transient  EM  analysis  the  TD-UTD 
solution  for  an  arbitrary  curved  conducting  wedge,  whose  development  was  indicated  last 
year  in  our  1994  JSEP  annual  report,  has  been  significantly  extended  to  arrive  at  a  TD-UTD 
slope  diffraction  solution  for  the  same  curved  wedge  geometry  when  it  is  excited  by  a  pulse 
with  a  rapid  spatial  variation  near  the  edge.  This  rapid  spatial  variation  can  occur,  for 
example,  when  the  pulsed  antenna  illuminating  a  wedge  exhibits  a  null  close  to  the  edge  of 
the  wedge.  The  transient  response  by  the  wedge  in  this  case  is  caused  by  the  large  slope 
(derivative)  of  the  rapid  spatial  variation  of  the  pulsed  incident  field,  and  the  ordinary  TD- 
UTD  solution  reported  last  year  would  produce  a  negligible  response  in  this  instance  thus 
requiring  the  need  for  the  slope  TD-UTD  curved  wedge  diffraction  solution.  In  addition, 
^  TD-UTD  solution  has  been  developed  for  the  shadow  side  response  of  a  smooth  convex 
conducting  surface  illuminated  by  a  pulsed  incident  field.  The  analytical  development  of  TD- 
UTD  for  convex  surfaces  is  found  to  be  even  more  involved  than  that  for  the  wedge,  although 
the  analytical  TD-UTD  solutions  obtained  for  these  geometries  is  still  relatively  simple  to 
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use.  Presently,  the  TD-UTD  for  the  lit  region  of  the  convex  surface  is  being  developed. 
The  TD-UTD  is  being  developed  separately  for  the  lit  and  shadow  sides  of  the  convex 
boundary  because  the  corresponding  FD-UTD  solution  which  is  analytically  being  inverted 
to  the  TD  is  also  obtained  separately  for  the  lit  and  shadow  sides  for  reasons  of  convenience 
of  analysis,  but  the  TD-UTD  lit  and  shadow  representations  should  be  continuous  at  the 
shadow  boundary  just  as  they  are  in  the  FD-UTD.  Also,  an  asymptotic  HF  (or  FD-UTD  like) 
has  been  developed  for  an  impedance  wedge  with  arbitrary  source  and  observer  locations, 
and  also  for  elongated  smooth,  convex  conducting  objects  where  the  source  and  observer 
can  again  be  located  arbitrarily.  The  latter  two  FD  solutions  are  crucial  to  the  development 
of  important  hybrid  combination  of  numerical-asymptotic  methods  for  the  analysis/design 
of  material  treated  antenna  configurations  for  optimizing  antenna  performance,  and  to  the 
development  of  radiation/scattering  from  aerospace  configurations,  respectively.  Some  more 
discussion  of  this  hybrid  connection  is  presented  later  in  the  section  on  Hybrid  Studies  under 
JSEP. 

Work  has  also  progressed  to  extend  the  reflection  and  diffraction  of  general,  astigmatic 
Gaussian  beams  (GB’s)  by  smooth  curved  boundaries  with  an  edge  in  which  the  edge  is  no 
longer  parallel  to  the  principal  planes  of  the  surface  where  the  axis  of  the  incident  astigmatic 
GB  strikes  the  surface  or  its  mathematical  extension  past  the  edge  beyond  which  the  actual 
surface  ceases  to  exist.  This  novel  GB  technique  is  expected  to  be  useful  especially  in 
predicting  EM  phenomena  for  focused  fields  where  even  ray  methods  become  cumbersome 
or  inapplicable.  In  particular,  such  a  GB  technique  is  expected  to  be  highly  useful  in  the 
analysis/synthesis  of  large  aperture  antennas  such  as  general  reflectors  (which  can  be  shaped 
or  array  fed,  etc.),  planar  phased  array  antennas  and  lens  antennas. 

2.  Time  Domain  UTD  (TD-UTD) 

The  TD-UTD  concept  provides  a  new  way  of  examining  transient  electromagnetic  (EM) 
radiation,  scattering  and  diffraction  phenomena  which  has  been  analyzed  in  the  past  for 
only  some  special  canonical  geometries,  or  for  some  more  complicated  but  relatively  small 
radiating  configurations  by  brute  force  numerical  methods  (which  depend  on  a  particular 
excitation  pulse  shape).  In  particular,  the  TD-UTD  extends  the  ray  concept  of  the  frequency 
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domain  (FD)  based  uniform  geometrical  theory  of  diffraction  (UTD)  to  the  time  domain 
(TD).  Such  a  TD-UTD,  when  it  is  sufficiently  developed  will  become  as  powerful  as  the  FD 
based  UTD;  it  would  be  able  to  predict  the  transient  response  of  pulse  excited  large  complex 
radiating  objects  in  the  form  of  an  analytical  solution  given  in  a  physically  appealing  manner 
in  terms  of  TD  rays.  Since  the  FD  based  UTD  is  an  asymptotic  high  frequency  technique,  its 
counterpart  in  the  TD  would  provide  a  TD-UTD  which  is  accurate  for  early  to  intermediate 
times  with  respect  to  the  time  of  arrival  of  each  of  the  various  wavefronts  associated  with 
the  different  ray  contributions.  As  is  well  known,  the  early  to  intermediate  time  response  is 
produced  largely  by  the  local  features  of  the  radiating  object,  whereas  the  late  time  response 
is  largely  dominated  by  the  global  effects  (including  resonances).  The  development  of  TD- 
UTD  is  motivated  by  the  need  to  analyze  transient  EM  phenomena  directly  in  the  TD.  Such 
an  analytical  prediction  tool  would  be  extremely  useful  for  assessing  the  effects  of  short  pulse 
radars,  lightning  or  nuclear  explosions,  etc.  on  communication  systems;  it  will  also  be  useful 
in  remote  sensing  and  target  identification  applications. 

Significant  developments  have  been  made  in  the  TD-UTD  area  since  its  initiation  on 
the  JSEP  program  about  two  years  ago.  The  developments  to  date  in  the  TD-UTD  have 
been  made  possible  largely  due  to  the  novel  use  of  analytic  time  transforms  (ATT),  since 
the  TD-UTD  is  being  developed  by  analytically  transforming  the  corresponding  FD  based 
UTD  solution  into  TD.  It  is  important  to  note  that  the  use  of  conventional  inverse  Fourier 
or  Laplace  transformations,  in  which  time  variable  (t) is  real,  generally  leads  to  difficulties  in 
obtaining  a  TD-UTD  from  its  FD  based  UTD;  these  difficulties  are  further  compounded  in 
the  case  when  the  TD-UTD  must  be  found  from  a  FD-UTD  in  which  the  incident,  reflected 
or  diffracted  rays  pass  through  their  respective  caustics.  In  the  latter  case,  conventional 
Fourier/Laplace  inversion  techniques  require  one  to  essentially  find  the  Hilbert  transform 
which  is  generally  extremely  difficult  to  obtain  analytically  in  closed  form.  By  making  the 
time  variable  ( t )  complex,  it  is  possible  to  use  the  ATT  and  thereby  make  the  inversion 
to  the  TD  case  resemble  the  direct  Laplace  transform  in  the  conventional  FD  except  that 
the  usual  roles  of  the  time  and  frequency  variables  are  now  switched.  Thus,  the  present 
TD-UTD  development  makes  it  possible  to  deal  with  concave  or  saddle  shaped  incident 
wavefronts  (i.e.,  other  than  just  convex  wavefronts)  with  great  ease.  Even  though  the  ATT 
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is  used,  the  actual  real  time  TD-UTD  is  obtained  by  taking  the  real  part  of  the  analytic 
time  solution  for  the  TD-UTD.  Furthermore,  the  use  of  the  ATT  has  also  led  to  another 
important  development  in  that  it  allows  the  TD-UTD  response  to  a  time  impulsive  general 
astigmatic  incident  ray  field  to  yield  in  closed  form  the  TD-UTD  response  to  a  finite  energy 
(finite  bandwidth)  incident  pulse  via  complex  t  convolution  (as  opposed  to  conventional  reed 
time  convolution)  by  expanding  the  incident  pulse  in  a  set  of  generally  few  special  basis 
functions  whose  ATT  can  also  be  found  in  closed  form;  these  basis  functions  are  simple 
exponentials  in  the  “frequency  domain”  or  FD.  Many  useful  transient  incident  waveforms  or 
pulse  shapes  can  be  synthesized  by  three  to  ten  such  basis  functions.  This  result  is  important 
because  it  allows  the  convolution  with  a  known  analytic  time  impulse  response  in  closed  form; 
the  efficiency  and  advantages  gained  by  obtaining  a  closed  form  TD-UTD  impulse  response 
would  otherwise  be  lost  by  a  numerical  convolution  which  is  now  avoided. 

During  the  previous  year’s  effort  in  which  the  TD-UTD  development  was  initiated,  use 
of  the  ATT  made  it  possible  to  obtain,  for  the  first  time,  a  TD-UTD  for  an  arbitrary 
curved  perfectly  conducting  wedge  excited  by  an  arbitrary  astigmatic  time  impulsive  ray 
field.  The  results  of  this  work  were  summarized  in  the  1994  JSEP  annual  report  for  last 
year’s  effort.  During  the  present  contract  period  the  TD-UTD  has  been  extended  to  two 
additional  and  important  situations;  namely  to  include  the  slope  diffraction  by  an  arbitrary 
curved  conducting  wedge,  and  to  the  diffraction  by  an  arbitrary  smooth  convex  conducting 
surface  in  the  shadow  region,  respectively. 

The  TD-UTD  for  slope  diffraction  by  a  wedge  was  obtained  by  inverting  two  different 
corresponding  FD  based  UTD  solutions  by  Hwang  and  Kouyoumjian,  and  by  Verrutipong 
and  Kouyoumjian,  respectively,  The  slope  diffracted  TD-UTD  is  required  when  the  incident 
transient  ray  field  varies  rapidly  in  space  near  the  point  of  diffraction  on  the  edge;  this  is 
important  for  many  antenna  applications  such  as  pulsed  apertures  (or  slots)  in  edge  struc¬ 
tures,  and  also  cases  in  which  the  antenna  pattern  illuminating  an  edge  exhibits  a  rapid 
spatial  variation,  as  for  example,  near  a  pattern  null.  In  this  instance,  the  slope  of  the  TD- 
UTD  for  the  curved  wedge  will  exhibit  a  discontinuity  near  the  incident /reflected  shadow 
boundaries  that  can  be  removed  by  the  addition  of  the  new  slope  TD-UTD  just  as  in  the 
FD  case.  The  ATT  was  used  again  and  it  was  this  analytic  time  signal  representation  that 
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made  it  possible  to  obtain  the  slope  TD-UTD  for  an  arbitrary  wedge.  The  development 
of  the  slope  TD-UTD  for  the  arbitrary  curved  wedge  was  slightly  more  involved  than  the 
development  of  the  TD-UTD  for  the  arbitrary  curved  wedge.  While  the  slope  TD-UTD  is 
also  obtained  for  a  general  curved  wedge,  the  numerical  results  based  on  this  TD-UTD  are 
illustrated  in  Figures  1  and  2  for  a  straight  wedge  because  an  exact  TD  analytical  solution 
exists  for  this  special  case  [1]  and  hence  it  can  be  used  as  a  reference  solution  to  validate 
the  more  general  TD-UTD  result.  The  field  incident  on  the  wedge  in  these  examples  is  time 
impulsive  spherical  wave  and  it  is  produced  by  an  electric  or  magnetic  current  moment  (with 
a  step  current  moment  behavior  in  time)  which  is  aimed  so  that  the  pattern  null  along  its 
axis  is  close  to  the  edge  to  produce  a  rapid  spatial  variation  of  the  incident  field  at  the  edge. 
It  is  seen  that  the  slope  TD-UTD  development  based  on  Hwang  and  Kouyoumjian  slope 
FD-UTD  (denoted  by  “Hwa”  in  the  figures)  behaves  slightly  better  for  intermediate  to  later 
times  with  respect  to  the  arrival  of  the  diffracted  wavefront,  while  the  other  slope  TD-UTD 
based  on  the  Verrutipong-Kouyoumjian  slope  FD-UTD  (denoted  by  “Ver”  in  the  figures)  is 
slightly  better  for  early  times.  Nevertheless,  both  TD-UTD  results  are  sufficiently  accurate 
in  general  in  the  early  to  intermediate  times  as  they  should  be. 

The  diffraction  by  a  smooth  convex  surface  is  far  more  involved  to  invert  into  the  time 
domain.  At  the  present  time,  a  significant  breakthrough  has  been  made  which  allows  ATT, 
or  analytic  time  signal  based  representation  to  provide  for  the  first  time  a  TD-UTD  for 
describing  the  fields  on  the  shadow  side  of  an  arbitrary  smooth  convex  perfectly  conducting 
boundary.  Work  is  in  progress  to  obtain  a  similar  TD-UTD  result  for  the  same  configuration 
on  the  lit  side  of  the  shadow  boundary;  this  lit  zone  analysis  appears  to  be  even  more 
challenging  than  the  shadow  side.  The  excitation  is  time  impulsive  arbitrary  astigmatic 
wavefront  that  is  incident  on  the  convex  surface.  It  is  found  for  the  shadow  side  that  the 
TD-UTD  version  of  the  FD-UTD  transition  function  can  be  represented  separately  for  the 
early  times  and  the  intermediate  (or  later)  times  with  respect  to  the  time  of  arrival  of 
the  diffracted  wavefront.  Furthermore,  these  two  representations  smoothly  overlap  thereby 
yielding  a  useful  TD-UTD  result.  It  turns  out  that  the  early  time  TD-UTD  representation 
comes  from  an  asymptotic  expansion  in  the  FD  requiring  a  large  number  of  so-called  creeping 
wave  mode  terms  to  be  summed  in  the  TD;  this  sum  has  been  converted  to  a  sum  on  a  small 
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Electric  Dipole  and  PEC  Wedge,  Diffraction 


Figure  1:  The  TD-UTD  impulse  response  is  compared  with  an  exact  result.  The  infinitesimal 
electric  current  has  a  unit  step  time  behavior  which  approximately  illuminates  the  wedge 
with  an  impulsive  spherical  wave. 


Electric  Field  Electric  Field 


Figure  2:  The  two  versions  of  the  TD-UTD  slope  diffraction  are  compared  with  an  exact 
result.  The  infinitesimal  magnetic  current  has  a  unit  step  time  behavior  which  approximately 
illuminates  the  wedge  with  an  impulsive  spherical  wave. 


number  of  terms  with  the  remainder  evaluated  in  closed  form  via  an  integral  approximation. 
For  later  times,  the  FD-UTD  transition  function  behaves  in  powers  of  frequency  yielding 
a  TD-UTD  representation  in  inverse  powers  of  time  (or  valid  for  intermediate  and  later 
times)  which  becomes  more  important  near  the  shadow  boundary.  Breaking  the  TD-UTD 
into  the  two  representations  which  turn  out  to  smoothly  overlap  in  the  shadow  region  yields 
a  tractable  TD-UTD  solution;  it  does  not  appear  possible  at  this  time  to  obtain  a  single 
representation.  Figures  3  and  4  illustrate  this  remarkable  overlap  of  the  two  representations 
for  the  TD-UTD  diffraction  by  a  circular  cylinder.  A  non-uniform  TD  solution  has  been 
obtained  previously  be  some  researchers  for  the  special  case  of  a  circular  cylinder;  in  contrast 
the  present  TD-UTD  is  a  uniform  solution  which  is  valid  for  an  arbitrary  convex  surface. 

Currently  work  is  in  progress,  as  mentioned  earlier,  to  complete  the  TD-UTD  develop¬ 
ment  for  the  diffraction  by  a  smooth  convex  surface  by  extending  it  into  the  lit  region.  Work 
has  also  been  initiated  to  develop  a  TD  version  of  the  equivalent  current  method  (ECM) 
as  well  as  of  the  physical  theory  of  diffraction  (PTD)  for  edged  bodies  to  augment  the  TD- 
UTD  edge  diffraction  solution  in  those  special  cases  where  there  is  an  overlap  of  ray  caustic 
and  shadow  boundary  regions  in  the  FD  so  that  the  FD-UTD  is  not  directly  applicable  and 
neither  is  the  corresponding  TD-UTD  for  that  reason. 

An  additional  related  study  which  was  undertaken  deals  with  the  development  of  an 
analytical  closed  form  prediction  of  early  time  transient  radiation  from  a  parabolic  reflector 
illuminated  by  a  pulsed  feed  antenna  at  the  focus.  The  feed  antenna  generates  a  pulsed 
spherical  wave  in  this  model.  It  is  seen  that  if  the  spherical  wave  from  the  feed  is  taken 
to  exhibit  a  step  function  behavior  in  time,  then  the  early  time  radiation  from  the  reflec¬ 
tor  exhibits  an  impulse  like  behavior  on  boresight.  Such  an  impulse  like  radiation  from  a 
parabolic  reflector  was  predicted  originally  by  C.  E.  Baum  and  his  co-workers  [2,  3]  when 
the  feed  was  a  TEM  four  wire  transmission  line  configuration;  they  were  able  to  predict  the 
behavior  of  the  fields  in  the  two  principal  planes  analytically  in  closed  form.  The  present 
analysis  proceeds  in  a  fashion  which  is  somewhat  different  than  the  one  given  previously  by 
Baum,  et  al.  [3]  and  it  provides,  to  within  the  model  chosen,  an  analytic  closed  form  solution 
for  the  fields  in  any  plane  (rather  than  only  the  principal  planes)  both  near  and  far  from  the 
reflector;  furthermore,  it  is  valid  for  a  more  general  spherical  wave  feed  illumination  and  it 
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Figure  4:  The  special  function  Fa  (H,t)  when  the  observer  is  in  the  shadow  region  and 
S  =  1.  The  polarization  is  soft.  Comparison  between  (a)  the  “early  time”  creeping  wave 
mode  series  and  (b)  the  “late  time”  representation. 
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can  be  modified  without  any  additional  complications  to  deal  with  an  offset  parabolic  reflec¬ 
tor.  In  this  sense,  the  present  work  extends  or  generalizes  the  earlier  work  on  this  topic.  A 
motivation  for  this  study  is  that  is  would  be  useful  in  the  design  of  pulsed  reflector  antennas 
which  can  radiate  concentrated  energy  for  use  in  an  impulse  like  radar  or  a  short  pulse  radar 
as  might  be  needed  in  target  identification  or  remote  sensing  application.  The  analytical 
closed  form  result  for  the  transient  radiation  by  a  parabolic  reflector  illuminated  by  a  step 
like  transient  spherical  wave  from  the  feed  can  be  differentiated  with  respect  to  time  to  pro¬ 
vide  the  transient  response  of  the  reflector  to  a  time  impulsive  spherical  wave  illumination 
from  the  feed;  in  the  latter  case,  the  reflector  radiates  a  doublet  like  field  on  boresight.  In 
addition,  the  transient  radiation  from  the  reflector  when  it  is  illuminated  by  the  feed  with 
a  transient  spherical  wave  and  a  given  pulse  shape  can  be  synthesized  via  convolution  with 
the  closed  form  impulse  response  (i.e.,  response  of  the  reflector  to  a  time  impulsive  spherical 
wave  from  the  feed).  An  example  of  the  impulse  like  radiation  from  the  parabolic  reflector 
when  it  is  illuminated  by  a  transient  step  (function)  type  spherical  wave  is  shown  in  the  more 
general  <j>  =  45°  plane  (rather  than  the  principal  planes  at  (f>  =  0°;  90°)  in  Figures  5  and  6, 
for  two  different  distances  from  the  reflector,  respectively.  On  boresight,  the  field  behaves 
more  and  more  like  an  impulse  as  one  moves  farther  away  from  the  reflector;  it  is  noted  that 
this  impulse  like  character  is  evidenced  by  the  fact  that  the  net  height  of  this  impulse  like 
field  on  boresight  is  independent  of  the  distance  away  from  the  reflector!  The  pulse  width  of 
the  impulse  like  transient  field  (on  boresight)  decreases  with  distance  but  its  pulse  amplitude 
increases  to  maintain  the  same  area  under  the  pulse;  however,  the  additional  range  effect 
(or  decrease  in  the  field  due  to  an  algebraic  range  factor)  exactly  nullifies  this  increase  in 
the  pulse  height  giving  a  net  field  amplitude  which  stays  constant  with  distance  from  the 
reflector  along  its  axis  (boresight).  Figures  5  and  6  also  indicate  the  manner  in  which  the 
impulse  like  behavior  is  lost  as  one  moves  off  boresight.  The  level  of  the  field  off  boresight  in 
the  <j>  =  45°  plane  is  seen  to  lie  somewhere  in  between  its  two  different  values  in  the  principal 
<t>  =  0°  and  4>  =  90°  planes,  respectively.  The  general  trend  of  the  field  behavior  in  the  two 
principal  planes  (<£  =  0°;  (f>  =  90°  planes)  is  seen  to  be  similar  to  that  in  the  <$>  =  45°  plane 
(illustrated  here  in  Figures  5  and  6),  or  for  that  matter  in  any  arbitrary  (f)  =  constant  plane. 
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3.  Gaussian  Beam  Methods 


Development  of  Gaussian  beam  (GB)  based  methods  were  initiated  on  our  JSEP  for  ana¬ 
lyzing  a  class  of  high  frequency  radiation/scattering  phenomena  that  cannot  be  analyzed 
efficiently,  or  in  a  tractable  fashion,  by  other  approaches  including  ray  methods.  In  particu¬ 
lar,  difficulties  often  arise  due  to  the  presence  of  caustics  in  the  ray  picture  when  employing 
ray  methods  for  describing  fields  in  regions  of  strong  focusing,  and  also  in  situations  where 
an  extremely  large  number  of  rays  are  required  to  describe  the  fields;  GB’s  can  efficiently 
handle  the  above  types  of  situations  because  they  yield  bounded  and  well  behaved  fields  in 
regions  of  focusing  where  ray  methods  become  singular.  With  this  view  in  mind,  a  GB  based 
expansion  was  developed  which  employs  general  astigmatic  GB’s  to  represent  the  fields.  The 
coefficients  of  this  expansion  can  be  found  approximately  but  accurately,  in  closed  form  in 
a  simple  manner,  or  quite  efficiently  by  a  numerical  Galerkin  procedure  if  desired.  This 
expansion  is  more  efficient  than  other  GB  expansions  which  employ  higher  order  Gauss- 
Hermite/Gauss-Laguerre  beam  modes;  also,  this  is  more  convenient  and  efficient  to  use  than 
a  Gabor  type  GB  expansion.  Such  a  GB  expansion  can  be  used  to  represent  an  incident  field, 
for  example,  and  the  radiation  or  scattering  phenomena  resulting  from  an  interaction  of  this 
incident  field  with  a  conducting  object  containing  an  edge  can  be  handled  by  using  essentially 
a  beam-optical  theory  developed  here.  In  particular,  such  a  theory  which  resembles  an  exten¬ 
sion  of  ray  methods  to  deal  with  general  astigmatic  GB  excitation  has  been  developed  which 
provides  closed  form  expressions  for  reflection  and  diffraction  of  general  astigmatic  GB’s  by 
a  perfectly  conducting  edge  in  an  otherwise  smooth  curved  surface.  This  development  was 
described  in  the  previous  1994  JSEP  annual  report.  However,  the  edge  diffraction  of  GB’s 
in  that  previous  work  was  restricted  to  edges  which  were  parallel  to  the  principal  directions 
of  a  local  paraboloidal  approximation  to  the  actual  surface  (or  its  analytic  extension)  at 
the  point  of  incidence  of  a  GB  near  the  actual  edge.  That  restriction  has  recently  been 
removed  such  that  the  edge  can  be  arbitrarily  located  but  the  edge  contour  must  lie  locally 
in  a  plane;  however,  it  appears  that  even  the  restriction  to  the  latter  case  of  a  planar  edge 
contour  can  be  removed  heuristically  and  is  under  investigation.  It  is  noted  that  most  edge 
contours  in  practical  antenna  applications  are  locally  planar.  The  extension  to  handle  the 
more  arbitrarily  located  planar  edge  was  developed  recently  in  a  tractable  fashion  by  using 
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asymptotic  evaluation  of  the  radiation  integral  for  the  scattered  field  when  the  illumination 
is  a  general  astigmatic  Gaussian  beam.  Such  an  asymptotic  evaluation  is  far  more  tractable 
for  this  more  complicated  edge  geometry  than  the  approach  used  in  the  earlier  development 
for  the  case  when  the  edge  was  parallel  to  the  local  surface  principal  directions  mentioned 
above;  however,  the  asymptotics  results  in  a  complex  stationary  point  which  must  be  found 
when  implementing  such  a  solution  in  analyzing  practical  phenomena  involving  GB  diffrac¬ 
tion  by  relatively  arbitrary  curved  edges.  It  appears  at  present  that  the  required  parameters 
in  the  solution  related  to  this  complex  stationary  point  can  be  found  approximately  but 
accurately  in  closed  form  in  which  case  the  implementation  of  this  solution  is  expected  to  be 
quite  efficient.  Work  is  continuing  to  study  this  situation.  Upon  the  completion  of  this  part 
of  the  study,  applications  of  this  solution  are  then  planned  to  be  initiated,  before  the  end  of 
this  year,  to  deal  in  turn  with  the  development  of  GB  synthesis  methods  for  electrically  large 
antennas.  In  particular,  the  fast  synthesis  of  reflector  antenna  feeds  and  shaped  reflectors, 
respectively,  will  be  investigated  to  provide  multiple,  spot  or  contour  beams  in  a  far  more 
efficient  and  physically  appealing  manner  than  is  currently  possible.  The  development  of  GB 
expansions  for  large  planar  phased  arrays  is  also  being  investigated.  Both,  the  use  of  planar 
phased  arrays,  and  the  use  of  large  shaped  reflectors  are  expected  to  play  a  very  significant 
role  in  the  future  for  both,  military  and  commercial  satellite  antenna  systems  in  wireless 
communication  applications.  Eventually,  it  is  also  proposed  to  develop  GB  techniques  for 
dealing  with  the  complex  phenomena  of  radiation/transmission  of  fields  of  large  antennas 
through  radomes  (including  frequency  selective  surfaces  (FSS)),  and  with  millimeter  wave 
systems.  It  is  noted  that  ray  techniques  are  not  directly  applicable  for  analyzing  many  of 
these  antenna  phenomena;  on  the  other  hand  physical  optics  or  integral  methods  are  found  to 
be  very  cumbersome  numerically.  The  GB  techniques  being  developed  here  are  not  expected 
to  suffer  these  limitations. 

4.  Diffraction  by  material  edges/junctions 

The  study  of  diffraction  by  material  bodies  (semi-transparent  or  opaque)  with  edges  and/or 
junctions  is  important  due  to  the  rapid  development  of  composite  materials.  The  analytical 
study  of  this  class  of  diffraction  phenomena  is  extremely  complex  and  so  far  only  a  very 


22 


limited  number  of  canonical  shapes  has  been  studied  rigorously  [4]- [9].  This  list  is  by  no 
means  complete  and  only  includes  few  papers  related  to  this  topic.  In  order  to  further  develop 
this  area  of  research,  it  is  necessary  to  use  some  combination  of  analytical  and  numerical 
techniques.  This  approach  allows  us  to  study  more  general  shapes  and  also  it  is  possible  to 
consider  fairly  arbitrary  material  properties.  The  approach  that  was  used  in  this  study  is  the 
Method  of  Moments/Green’s  function  scheme  [10],  where  the  Green’s  function  of  a  canonical 
geometry  is  obtained  analytically  (if  possible)  and  any  perturbation  of  the  canonical  geometry 
is  handled  by  the  Method  of  Moments.  The  canonical  geometry  considered  in  this  study  is  a 
material  cylinder  of  arbitrary  shape  and  material  properties  which  is  located  near  or  attached 
to  a  wedge.  Note  that  the  faces  of  the  wedge  satisfy  impedance  boundary  conditions,  which 
provides  the  means  for  studying  objects  whose  surfaces  are  non-perfectly  conducting,  such 
as  coated  conductors,  conductors  with  surface  roughness  or  corrugations.  The  very  general 
configuration  of  a  wedge  was  chosen  for  its  versatility  in  representing  a  wide  variety  of 
geometries.  Most  of  the  details  of  this  work  are  discussed  in  the  the  Hybrid  Studies  unit. 
It  is  important  to  note  that  the  work  that  is  relevant  to  this  diffraction  studies  unit  has  to 
do  with  the  development  of  the  two-dimensional  Green’s  function  for  the  impedance  wedge 
shown  in  Figure  7.  It  turns  out  that  the  key  step  in  the  present  hybrid  approach  is  to  obtain  a 
Green’s  function  for  the  impedance  wedge  that  is  numerically  efficient.  The  Green’s  function 
thus  developed  also  provides  a  useful  physical  insight  into  the  various  scattering  mechanisms 
that  exist  because  it  allows  the  individual  components  of  the  fields  to  be  isolated;  such  as 
field  components  that  diffract  from  the  tip,  reflect  from  the  faces,  and  propagate  along  the 
faces,  etc. 

The  Green’s  function  obtained  in  this  study  for  the  geometry  in  Figure  7  is  valid  for  all 
source  and  observation  locations.  There  have  been  some  developments  for  the  wedge  Green’s 
function  in  the  past,  all  of  which  have  had  certain  restrictions  placed  on  them  [8]— [9] .  The 
development  of  the  present  Green’s  function  for  the  impedance  wedge  depicted  in  Figure  7 
is  based  on  a  method  originally  developed  by  Clemmow  [11].  The  keys  steps  in  this  analysis 
are  (1)  obtain  the  plane  wave  response  of  the  impedance  wedge,  (2)  use  analytic  continuation 
to  extend  the  results  of  the  plane  wave  response  for  complex  angles  of  incidence,  (3)  expand 
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the  incident  field  from  the  line  source  into  a  spectrum  of  plane  waves  and  (4)  represent  the 
line  source  response  of  the  wedge  as  a  spectrum  of  plane  waves. 

The  line  source  response  for  the  wedge  can  be  written  as 


U{p,Po)  =  ^~  f  Up(<t>,a)eikpoCos^°-a)da. 
47r  J 


where  S( 7t/2)  is  an  Steepest  Descent  Path  (SDP)  through  a  =  ir/2. 

The  total  field  can  be  represented  as  the  Sommerfeld  integral 

Up(<p,  =  (  U{ w  +  mr/2  -  <p)e~ikpcosudu>,  (2) 

Z7TZ  J 
7 

where  7  is  the  twofold  Sommerfeld  contour.  The  term  U  is  given  by 

fT,  ,  V{u>  +  mr/2-<P)  sm{<p'/n) 

U{u  +  mr2-(p)  =  — — - — - — - — - ~rrz.  (3) 

nW{mr/2  —  <p')  cosf5^)  —  cos(^-) 

where  $  represents  the  product  of  four  Maliuzhinets  functions  ip  with  different  arguments. 
The  resulting  expressions  are,  in  general,  extremely  complicated  and  computationally  inten¬ 
sive  to  evaluate.  It  is  therefore  important  that  every  effort  be  made  to  make  them  as  efficient 
as  possible  to  evaluate  numerically.  For  certain  special  cases,  such  as  when  the  observation 
or  source  points  are  either  extremely  far  or  extremely  close  to  the  wedge  vertex,  it  will  be 
possible  to  asymptotically  approximate  these  integrals  with  a  closed  form  expression.  This 
will,  however,  not  be  the  case  in  general  and  the  integrals  have  to  be  evaluated  numerically. 
To  numerically  evaluate  such  integrals,  the  contour  of  integration  is  deformed  to  one  for 
which  the  integrand  decays  very  rapidly  as  the  magnitude  of  the  variable  of  integration  in¬ 
creases.  This  allows  the  range  of  integration  to  be  truncated.  One  of  two  types  of  contours 
of  integration  are  chosen,  either  an  SDP  or  a  Sommerfeld  contour.  When  the  distance  pa¬ 
rameter  is  larger  than  a  certain  specified  value,  the  SDP  is  used,  otherwise  the  Sommerfeld 
contour  is  used.  In  addition  to  deforming  the  contours  of  integration,  the  singularities  of  the 
integrands  axe  explicitly  separated  and  evaluated  in  closed  form. 

The  results  for  the  development  of  the  Green’s  function  were  verified  numerically  by 
comparing  it  with  a  MM  reference  solution.  The  reference  solution  was  developed  by  for¬ 
mulating  an  integral  equation,  which  utilizes  a  ffee-space  Green’s  function,  to  determine  the 
scattering  from  an  impedance  wedge.  The  integral  equation  was  then  solved  using  a  MM 
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technique  that  utilized  pulse  basis  functions  and  point  matching.  If  either  the  source  or 
observation  point  is  less  than  0.25A  from  the  vertex,  the  form  of  the  Green’s  function  which 
was  developed  for  the  case  when  one  of  the  distance  parameters  is  small  is  used,  otherwise 
the  form  developed  for  the  case  when  it  is  large  is  used.  For  the  MM  solution,  50  pulse  basis 
functions,  each  of  width  0.1  A,  is  placed  on  each  of  the  wedge  faces.  Since  we  are  modeling  an 
infinite  structure,  enforcing  the  impedance  boundary  condition  only  along  5  A  long  sections 
of  the  wedge  faces  could  be  a  problem  because  it  creates  spurious  endpoints  that  could  be 
points  of  diffraction  for  any  unattenuated  surface  waves  propagating  away  from  the  vertex. 
Therefore,  only  for  the  sake  of  comparison  with  the  results  based  on  the  MM  solution,  we 
will  be  considering  impedance  surfaces  that  are  lossy  so  that  any  surface  waves  that  are 
excited  will  be  sufficiently  attenuated  when  they  propagate  along  the  wedge  face  and  there 
will  be  a  negligible  interaction  with  these  endpoints. 

First,  we  will  consider  the  case  when  the  line  source  is  at  pa  =  0.2  A  and  the  observation 
distance  is  p  =  0.75  A;  therefore  we  are  using  the  form  of  the  Green’s  function  developed  for 
the  small  distance  parameter.  Two  different  wedge  configurations  are  considered  where  the 
exterior  wedge  angle  may  be  330°  or  210°.  The  results  for  the  two  wedges  with  the  surface 
impedance  values  Zx  =  37.7  —  z566fl  /sq.  and  Z2  =  189  —  z377Q  /sq.  are  shown  in  Figures  8 
and  9  for  both  TM2  and  TE2  polarizations.  These  are  bistatic  plots  of  the  total  field  vs.  the 
observation  angle  <f>. 

In  all  of  these  cases  the  agreement  between  the  Green’s  function  and  MM  solution  is 
excellent.  When  the  imaginary  part  of  the  surface  impedance  is  negative,  that  surface  can 
support  surface  waves  for  a  TE2  field  but  not  for  a  TM2  field,  whereas  when  the  imaginary 
part  of  the  surface  impedance  is  positive,  it  is  the  TM2  field  that  can  excite  surface  waves 
and  not  the  TE2  field.  The  presence  of  surface  waves,  when  they  are  present,  is  clearly  seen 
by  the  sharp  increase  of  the  fields  as  the  observation  point  approaches  the  wedge  faces,  with 
the  stronger  surface  wave  being  along  face  1.  This  is  to  be  expected  since  the  magnitude  of 
the  imaginary  part  of  Zx  is  larger  than  that  for  Z2  and  face  1  has  less  loss.  Also,  note  that 
all  of  the  plots  for  the  Green’s  function  are  continuous  across  the  many  shadow  boundaries 
that  exist  for  these  solutions. 
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Figure  8:  Bistatic  plot  comparing  the  Green’s  function  solution  to  a  moment  method  solution. 
The  parameters  are  pa  =  .2A,  p  =  .75A,  <f>0  =  120°,  Z\  =  37.7  —  i566  fi/sq.,  Zi  =  189  — 
i377  fi/sq.,  and  the  exterior  wedge  angle  is  330°. 
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GREEN'S  FUNCTION 


Figure  9:  Bistatic  plot  comparing  the  Green’s  function  solution  to  a  moment  method  solution 
The  parameters  are  pQ  =  .2A,  p  =  .75A,  4>0  -  120°,  Zl  =  37.7  -  z‘566  Q/sq.,  Z2  =  189  - 
i377  fl/sq.,  and  the  exterior  wedge  angle  is  210°. 


RADIATING 

PATCH 


Figure  10:  Two-dimensional  configuration  for  a  microstrip  patch  antenna  on  a  semi-infinite, 
grounded,  dielectric  substrate  (a),  and  the  equivalent  models  for  the  TE2  case  (b),  and  the 
TM2  case  (c). 

The  two-dimensional  (2-D)  wedge  Green’s  function  is  also  applicable  for  studying  edge 
effects  on  the  radiation  pattern  of  antennas.  One  of  the  types  of  antennas  that  can  be 
effectively  modelled  with  this  Green’s  function  is  a  microstrip  patch  mounted  on  a  grounded 
dielectric  substrate,  an  example  of  which  is  shown  in  Figure  10(a). 

Here,  the  substrate  is  terminated  at  one  end  leaving  an  edge  that  can  give  rise  to  a 
point  of  diffraction.  The  fields  that  diffract  from  that  edge  can  have  a  significant  impact  on 
the  radiation  pattern  and,  hence,  degrade  the  antenna’s  performance.  The  wedge  Green’s 
function  allows  us  to  isolate  the  effect  of  a  single  edge  on  the  radiation  pattern  and  it  also 
provides  the  means  to  study  how  the  different  field  components  interact  with  the  edge. 
Although  our  model  is  2-D,  it  can  be  effectively  used  to  analyze  a  3-D  configuration  if  we 
use  the  TM2  and  TE2  polarizations  in  the  2-D  case  to  analyze  the  H-plane  and  E-plane 
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patterns,  respectively,  for  the  3-D  case.  The  configuration  in  Figure  10(a)  is  modelled  by  an 
impedance  wedge  with  an  exterior  wedge  angle  of  2tt,  the  grounded  dielectric  substrate  is 
replaced  by  the  equivalent  impedance  surface  Zs  for  the  top  face  of  the  wedge,  the  bottom 
face  remains  a  PEC,  and  the  radiating  patch  is  replaced  by  its  equivalent  currents;  either 
two  magnetic  line  sources  for  the  TE2  case  (Figure  10(b)),  or  a  sheet  of  electric  current  for 
the  TM2  case  (Figure  10(c)). 

One  of  the  advantages  of  using  the  wedge  Green’s  function  in  this  analysis  is  that  the 
different  field  components  can  be  explicitly  separated  and  studied  individually.  This  can 
be  a  valuable  tool  in  providing  physical  insight  into  the  various  mechanisms  involved  in 
the  radiation  and  scattering  from  structures  and  hence  into  their  control  to  achieve  desired 
performance.  In  determining  the  radiation  pattern  of  the  patch  antenna  in  the  presence  of 
the  edge,  there  are  three  field  components  that  should  be  considered.  The  direct  space  wave 
represents  the  field  that  is  radiated  by  the  patch  and  travels  directly  to  the  observation  point 
without  interacting  with  the  edge,  and  then  there  is  a  component  of  this  space  wave  that 
diffracts  from  the  edge.  Since  there  is  an  impedance  surface,  surface  waves  can  be  excited 
which  travel  away  from  the  patch  along  the  surface  and  the  one  that  travels  toward  edge 
can  diffract  and  contribute  to  the  total  field  seen  at  the  observation  point.  For  the  structure 
shown  in  Figure  10  with  Zs  —  —  i34fi/sq.  and  d  =  0.184A,  the  E-plane  radiation  patterns 
are  shown  in  Figure  11  for  the  direct,  diffracted  space  wave  (DIFF),  diffracted  surface  wave 
(SW),  and  total  fields.  The  direct  field  is  the  field  that  exists  for  the  antenna  mounted 
on  an  infinite  substrate,  no  edge  being  present.  The  effect  of  the  edge  is  to  give  rise  to 
a  point  of  diffraction  that  causes  a  considerable  amount  of  radiation  in  the  back  region, 
180°  <  4>  <  360°.  The  main  contribution  to  the  scattered  fields  in  this  region  comes  from 
the  space  wave  that  diffracts  from  the  edge,  whereas  the  surface  wave  contribution  is  almost 
negligible.  However,  if  the  surface  impedance  were  higher,  as  a  result  of  a  thicker  dielectric 
which  is  modeled  by  this  surface  impedance,  for  example,  then  the  contribution  from  the 
surface  wave  would  be  more  noticeable.  Note  that  the  Green’s  function  presented  here 
can  be  made  more  general  by  replacing  the  Leontovich  boundary  conditions  by  Generalized 
Impedance  boundary  conditions  [12]. 
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Figure  11:  TE2  radiation  pattern  for  the  configuration  shown  in  Figure  10  with  Zs  = 
— i34fl/sq.,  d  =  0.184A,  and  L  =  0.3A.  The  individual  field  terms  are;  the  total  field  in 
the  presence  of  the  edge  (TOTAL),  direct  field  term  that  does  not  include  the  edge  (DI¬ 
RECT),  diffracted  space  wave  (DIFF),  and  diffracted  surface  wave  (SW). 
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5.  Extensions  of  UTD  for  Conducting  Surfaces 


A  UTD  solution  has  been  developed  for  predicting  the  near  fields  (but  not  extremely  near 
fields)  and  far  fields  of  antennas  that  are  placed  close  to  a  smooth  perfectly-conducting  convex 
surface;  these  results  reduce  to  the  ones  previously  obtained  [13]  for  antennas  located  directly 
on  the  smooth,  perfectly-conducting  convex  surface.  Work  is  also  in  progress  to  obtain  an 
extended  UTD  solution  for  source  (antenna)  and  observer  (field  point)  which  are  both  close 
to  a  smooth  perfectly-conducting  convex  surface;  this  solution  is  required  for  numerically 
blending  into  the  other  solution  valid  for  an  antenna  close  to  (but  for  an  observer  not  close  to) 
the  smooth  convex  boundary,  which  has  recently  been  developed  as  mentioned  above.  Also, 
this  extended  UTD  solution,  which  is  currently  under  development,  is  expected  to  reduce 
to  the  UTD  solution  obtained  previously  [14]  for  source  and  observer  that  are  both  directly 
on  a  smooth  convex  boundary.  The  reason  for  obtaining  these  different  UTD  solutions  (or 
different  UTD  representations)  for  different  regions  of  source  and  observer  location  is  that 
such  different  representations  can  be  obtained  in  terms  of  UTD  transition  functions  that 
are  dependent  only  on  a  single  parameter,  and  which  are  therefore  also  well  tabulated.  On 
the  other  hand,  a  general  solution  valid  for  arbitrary  source  and  observer  location  involves 
the  use  of  special  UTD  transition  functions  which  are  dependent  on  several  parameters,  and 
therefore  cannot  be  easily  tabulated  as  are  universal  functions  of  a  single  parameter.  It  is  for 
this  reason  that  a  general  UTD  solution  valid  for  arbitrary  source  and  observer  location  is 
not  being  employed;  instead,  different  representations  of  such  a  general  solution,  which  are 
simpler  to  use  from  a  numerical  computation  point  of  view  and  which  are  to  remain  valid  for 
different  range  of  source  and  observer  locations  are  sought  in  this  study.  It  is  expected  that 
the  extended  UTD  solution  for  source  and  observer  both  very  close  to  the  convex  surface 
will  be  completed  by  the  end  of  this  year.  Also,  work  has  progressed  significantly  to  obtain 
an  asymptotic  high  frequency  solution,  for  the  scattering  by  an  elongated  smooth  convex 
surface  with  an  endcap,  which  remains  valid  in  the  paraxial  region  of  the  elongated  surface, 
to  simulate  the  scattering  from  the  fuselage  of  an  aircraft /missile  in  which  the  end  cap 
models  the  bulkhead  for  the  nose  antenna  inside  a  radome  (which  is  assumed  transparent). 
This  asymptotic  solution  is  also  expected  to  be  completed  by  the  end  of  the  year.  The  main 
reasons  for  the  development  of  all  the  UTD,  or  asymptotic  solutions  described  under  this 
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subtopic  are  that  they  will  be  crucial  for  not  only  predicting  additional  important  effects  of 
aircraft /missile  structure  on  the  performance  of  antennas  located  on  such  complex  realistic 
structures  which  cannot  be  handled  presently,  but  also  be  crucial  for  predicting  radiation 
and  scattering  of  EM  waves  by  such  structures  via  a  hybrid  MM-UTD  approach  discussed 
in  slightly  more  detail  within  the  Hybrid  Studies  unit  of  JSEP. 

6.  Development  of  UTD  Solutions  for  Arbitrary  Convex  Con¬ 
ducting  Surfaces  with  Uniform  Material  Coating 

Recently,  highly  useful  UTD  solutions  have  been  developed  to  analyze  the  radiation  from, 
and  mutual  coupling  between  antennas  located  on  or  in  material  coated  smooth  perfectly 
conducting  convex  surfaces.  Actually,  the  antenna  element  chosen  in  this  UTD  development 
is  a  point  electric  or  magnetic  current  source;  however,  the  UTD  solutions  for  a  point  source 
excitation  can  be  used  to  directly  obtain  the  corresponding  UTD  solutions  for  an  arbitrary 
antenna  source  distribution  via  superposition.  Thus,  UTD  Green’s  functions  for  exterior 
radiation  and  surface  fields  have  been  developed  for  point  source  excitation  of  a  coated 
conducting  convex  surface.  It  is  assumed  that  the  material  coating  is  homogeneous  and 
of  uniform  thickness;  also  the  coating  is  assumed  to  be  thick  with  respect  to  the  radii  of 
curvature  of  the  convex  surface  which  in  turn  are  assumed  to  be  electrically  large.  These 
UTD  solutions  properly  reduce  to  those  obtained  previously  [13,  14]  also  under  JSEP  for  a 
perfectly-conducting  smooth  convex  surface,  in  the  limit  when  the  thickness  of  the  coating 
is  allowed  to  vanish.  Presently,  the  special  UTD  transition  functions,  which  involve  modified 
Fock  functions  for  this  coated  case,  are  being  coded  so  that  they  can  be  computed  very 
efficiently  making  these  UTD  solutions  not  only  physically  appealing  but  also  useful  in 
analyzing  the  performance  of  antennas  which  truly  conform  to  an  arbitrary  3-D  smooth 
perfectly  conducting  convex  surface  with  a  uniform  material  coating.  It  is  expected  that  the 
computations  of  the  special  UTD  transition  functions,  which  occur  in  these  UTD  Green’s 
functions  for  the  coated  convex  surface,  will  be  completed  before  the  end  of  this  year.  As 
such,  these  UTD  solutions  would  provide  for  the  first  time,  essentially  closed  form  analytical 
results  for  predicting  the  effects  of  surface  curvature  and  material  coating  on  the  performance 
of  antennas  located  on  such  curved  coated  boundaries.  Surface  ray  torsion  explicitly  shows  in 
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these  UTD  solutions  how  the  polarization  properties  of  antennas  on  or  in  coated  conducting 
surfaces  change  with  the  direction  of  radiation,  and  also  how  that  affects  antenna  mutual 
coupling.  If  the  antenna  source  distribution  is  not  known  exactly  or  even  approximately,  then 
these  UTD  Green’s  functions  can  be  employed  as  the  kernel  of  an  integral  equation  for  the 
unknown  antenna  source  distribution.  The  unknowns,  because  of  the  special  choice  of  the 
UTD  Green’s  functions,  reside  only  over  the  relatively  small  regions  of  the  antenna  source 
distribution  and  can  thus  be  found  very  efficiently  via  the  numerical  moment  method  (MM) 
thereby  leading  to  a  hybrid  MM-UTD  procedure.  Such  a  hybrid  scheme  is  particularly  useful 
for  the  analysis/design  of  large  but  obviously  finite  antenna  arrays  with  relatively  arbitrary 
antenna  elements  that  conform  to  an  arbitrarily  curved  3-D  smooth  perfectly  conducting 
surface  with  a  thin,  uniform  material  coating;  this  aspect  is  discussed  later  in  the  section  on 
Hybrid  Studies  under  JSEP. 
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I.  Tekin,  Graduate  Research  Associate  (Phone:  614/292-7981) 

D.  Torrungrueng,  Graduate  Research  Associate  (Phone:  614/292-7981) 

1.  Overview 

This  section  will  summarize  our  work  in  integral  equation  studies  from  September  1994  to 
September  1995.  This  year  represents  a  change  in  direction  for  the  Integral  Equation  Studies 
area.  In  particular,  over  the  past  few  years  our  research  has  concentrated  on  the  integral 
equation  and  method  of  moments  (MM)  analysis  of  artificial  media  [l]-[6]  and  extremely 
low  frequency  electromagnetic  shielding  [7].  This  work  is  now  over,  and  our  present  research 
involves  two  topics: 

1.  The  Time  Division  Method  of  Moments  (TDMM)  analysis  of  the  scattering  by  electri¬ 
cally  large  bodies. 

2.  Two  important  junction  problems  necessary  for  the  MM  analysis  of  arbitrary  bodies. 
Our  work  over  the  past  year  on  these  two  topics  is  briefly  summarized  below. 

2.  The  TDMM  Analysis  of  Electrically  Large  Bodies 

a.  Introduction 

In  the  MM,  the  current  on  a  body  is  approximated  by  an  N  term  expansion  of  the  form 

j*j"  =  f;jnFn  (4) 

n= 1 

where  the  Fn  are  a  series  of  N  known  basis  or  expansion  functions  and  the  J„  are  a  series 
of  N  unknown  coefficients.  By  defining  suitable  weighting  functions  and  applying  the  MM 
procedure,  the  /„  are  found  as  a  solution  of  the  order  N  matrix  equation 

[Z)I  =  V  (5) 
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where  \Z\  is  the  the  order  N  impedance  matrix,  V  is  the  length  N  voltage  vector,  and  /  is 
the  length  N  current  vector  which  holds  the  unknown  coefficients,  J„,  from  Equation  4. 

The  number  of  terms,  N,  which  must  be  kept  in  the  summation  of  Equation  4,  in  order 
that  JA  be  a  reasonable  approximation  to  the  true  current  J,  is  proportional  to  the  electrical 
size  of  the  body.  For  large  N,  the  major  computational  resources  necessary  to  carry  out  the 
MM  solution  are 

•  the  0(N2)  CPU  time  to  compute  the  N2  elements  of  the  [Z]  matrix. 

•  the  0(N2)  (complex)  memory  locations  to  store  the  N 2  elements  of  the  [Z]  matrix. 

•  the  0(N3)  CPU  time  to  solve  Equation  5  for  I  by  direct  methods  such  as  LU  decom¬ 
position. 

The  major  obstacle  in  applying  to  the  MM  to  electrically  large  bodies  is  the  0(1V3)  solution 
time.  Over  the  past  few  years  several  different  methods  have  been  developed  in  order  to 
reduce  this  CPU  time  including  iterative  solutions  of  the  matrix  equation,  the  use  of  the  fast 
multipole  method  (FMM)  to  reduce  matrix  storage  and  the  CPU  time  per  iteration,  spatial 
decomposition  techniques  in  which  a  large  body  is  split  into  smaller  parts,  combining  the 
MM  with  an  asymptotic  technique  such  as  the  GTD,  the  use  of  wavelet  and  other  special 
basis  function  to  obtain  a  sparse  matrix,  and  many  others  [8]. 

The  method  presented  here  is  termed  the  Time  Division  Method  of  Moments  (TDMM). 
Although  it  is  a  frequency  domain  solution  for  the  time  harmonic  current  on  a  body,  it  is 
helpful  to  visualize  it  as  the  current  induced  by  a  time  unit  step  plane  wave  excitation  of 
the  scatterer.  The  scatterer  will  be  analyzed  in  a  sequence  of  time  intervals  which  directly 
correspond  to  spatial  sections  on  the  body.  As  one  proceeds  for  one  time  interval  to  the 
next,  the  scatterer  current  can  only  change  on  a  single  section  of  the  body,  thus  allowing  us 
to  analyze  the  large  order  N  problem  as  a  sequence  of  smaller  order  N/P  problems,  where 
P  is  the  number  of  sections  into  which  the  scatterer  has  been  split. 

For  simplicity  the  method  will  be  described  for  the  2D  problem  of  TE  to  z  scattering  by 
a  perfect  electric  conducting  (PEC)  strip,  however,  it  should  be  applicable  to  any  smooth 
convex  body.  The  problem  is  to  find  the  time  harmonic  current  induced  on  the  strip  by  the 
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plane  wave 


H'z(x,y,t)  =  cJ*o(*co8*+».in  <t>i)ejUt  (6) 

incident  with  angle  <f>i  on  a  (PEC)  strip  of  width  W.  Here  ko  is  the  free  space  wavenumber 
and  u>  =  27 rf  is  the  radian  frequency. 

As  illustrated  in  Figure  12,  in  the  TDMM  method  it  is  helpful  to  think  of  the  strip  excited 
by  the  unit  step  plane  wave 


H'zu{x,y,t )  =  H'z(x,y,t)u[c(x  cos 4>,  +  y  sin  </>,)  +  t] 


(7) 


where  c  is  the  speed  of  light,  t  is  time,  and  the  unit  step  function 


u(x) 


0  x  <  0 
1  x  >  0 


(8) 


For  the  purposes  of  a  standard  MM  solution,  the  strip  is  segmented  into  N  sections  of  width 
d  =  W/N,  corresponding  to  the  N  expansion  functions  of  Equation  4.  For  the  purposes  of 
the  TDMM  solution  the  MM  basis  functions  are  grouped  into  P  sections  of  width  D  =  W/ P, 
and  with  Np  =  N/P  basis  functions  per  section.  The  advantage  of  the  TDMM  method  it 
that  it  is  only  necessary  to  solve  equations  of  order  N/P ,  as  opposed  to  N  for  a  standard 
MM  solution. 


b.  The  First  Sweep 

The  TDMM  method  begins  with  the  k  =  1  or  first  sweep  which  computes  the  current  caused 
by  the  unit  step  incident  plane  wave  as  it  first  washes  across  the  strip,  i.e.,  in  the  time 
interval 

0  <  t  <  AT  AT  =  W |  cos  (j>i\/c  (9) 

from  when  the  wavefront  first  hits  the  left  edge  of  the  strip  to  just  before  it  hits  the  right 
edge.  The  first  sweep  begins  by  computing  the  forced  response  current  caused  by  the  unit 
step  incident  plane  frozen  in  time  at 

t  =  6T~  6T  =  AT/P  (10) 

just  before  the  wavefront  hits  the  right  edge  of  section  p  =  1.  At  this  time  there  is  no  current 
on  sections  p  =  2, 3, . . . ,  P  of  the  strip,  and  thus  these  sections  can  be  replaced  by  free  space. 
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Figure  12:  For  the  TDMM  solution  the  PEC  strip  is  split  into  P  sections  of  N/P  unknowns 
per  section. 

As  illustrated  in  Figure  13(a),  the  current  on  sweep  1  section  1  is  that  induced  by  the  time 
harmonic  incident  plane  wave  H'z,  except  for  the  undesired  presence  of  the  right  edge  of 
section  1.  To  electromagnetically  remove  this  edge,  a  tapered  resistance  or  R-card  is  added 
to  the  right  edge  of  section  1.  The  R-card  has  a  resistance  which  tapers  smoothly  from  a 
small  resistance  at  its  left  edge  (junction  with  PEC  section)  to  a  large  resistance  at  its  right 
edge  (junction  with  free  space).  Typically,  the  R-card  width  will  be  about  a  wavelength,  but 
for  practical  reasons,  we  always  choose  the  R-card  to  be  an  integer,  Pr,  number  of  sections. 
In  Figure  13  the  R-card  is  shown  as  Pr  =  2  sections.  Omitting  the  details,  assume  that  the 
current,  J^J,  on  section  p  =  1  for  sweep  k  =  1  has  been  computed. 

The  next  step  for  sweep  1  is  to  compute  the  currents  on  sections  p  =  2,3,  ...,P.  As 
illustrated  in  Figure  13(b),  the  sweep  1  currents  on  section  p  >  1  are  produced  by  the 
superposition  of: 

•  the  time  harmonic  incident  field  H'(x,y,t) 

•  the  time  harmonic  currents,  J*,  previously  computed  on  sections 

1, 2, . . . ,  p  —  1 
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(b)  Sweep  1  Current  on  Section  p 

Figure  13:  Computation  of  the  sweep  1  currents  for  a  PEC  strip. 
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radiating  in  the  presence  of  PEC  section  p  plus  Pr  =  2  sections  of  tapered  R-card  to  the 
right  of  section  p.  Note  that  in  the  TDMM  method,  all  computations  are  frequency  domain 
or  time  harmonic.  The  concept  of  time  stepping  is  for  conceptual  purposes  only. 

To  illustrate  the  sweep  1  currents,  the  insert  in  Figure  14  shows  a  wave  incident  from  &  = 
135°  on  a  strip  of  width  W  —  12A.  The  standard  MM  solution  follows  that  of  Richmond  [9], 
and  employs  a  segment  size  of  d  =  0.1A,  thus  resulting  in  N  =  120  unknowns.  For  the  TDMM 
solution  the  basis  functions  are  grouped  into  P  =  10  sections  of  width  D  =  W/P  =  1.2A, 
and  containing  Np  —  N/P  =  12  unknowns  each  *.  The  R-card  occupies  Pr  =  2  sections, 
and  thus  is  of  width  Dr  =  PrD  —  2.4A.  Its  resistance  tapers  exponentially  from  0  to 
13,000  f 2/D.  Figure  14(a)  compares  the  standard  MM  (solid  line)  and  sweep  1  (dashed  line) 
solutions  for  the  magnitude  of  the  MM  current  vector,  |/„|,  n  =  1,2, ... ,  N.  For  comparison, 
the  figure  also  shows  the  physical  optics  current  (short  dashed  line).  Note  that  the  sweep 
1  currents  are  significantly  better  than  the  physical  optics  current,  since  they  include  the 
effects  of  diffraction  at  the  leading  edge  of  the  strip.  The  higher  frequency  oscillations  seen 
in  the  standard  MM  solution  are  not  seen  in  the  sweep  1  TDMM  solution.  This  is  not  an 
error  in  the  sweep  1  currents,  but  rather  is  a  result  of  the  fact  that  the  higher  frequency 
oscillations  are  caused  by  multiple  reflections  of  the  current  which  occur  at  later  times.  Since 
for  sweep  1  the  unit  step  incident  wave  has  not  yet  hit  the  right  or  trailing  edge,  the  sweep  1 
currents  should  be  precisely  those  of  a  PEC  half-plane  x  >  0.  Figure  14(a)  shows  the  exact 
half-plane  currents  (dotted  line)  which  are  virtually  identical  to  the  sweep  1  currents. 

c.  The  Second  and  Subsequent  Sweeps 

When  the  wavefront  hits  the  right  or  trailing  edge  of  the  strip,  it  will  produce  a  reflected 
current  which  will  modify  the  current  on  section  P,  then  section  P-1,  . . .,  and  finally  section 
1.  To  account  for  this  progression,  sweep  2  is  done  in  reverse  order.  As  illustrated  in  Figure 
15,  the  sweep  2  currents  on  section  p,  j£=2,  are  produced  by  the  superposition  of: 

•  the  time  harmonic  incident  field  H'(x,y,t) 

•  the  time  harmonic  sweep  1  currents,  J*,  previously  computed  on  sections  1, 2, . . .  ,p  —  1 

1  Actually  due  to  the  use  of  the  overlapping  piecewise  sinusoidal  basis  functions,  there  are  IV  =  119 
unknowns,  with  only  11  unknowns  in  the  last  section. 
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Figure  14:  A  comparison  of  the  standard  MM  and  the  sweep  k  =  1,2,10  TDMM  currents  on 
a  PEC  strip. 
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•  the  time  harmonic  sweep  2  currents,  previously  computed  on  sections  p  +  l,p  + 
2,...,P 

radiating  in  the  presence  of  PEC  section  p.  Note  that  no  R-cards  are  required  on  the  the 
second  or  subsequent  sweeps  since  previously  computed  currents  enforce  continuity  of  current 
at  both  edges  of  section  p  (except  for  section  1  and  P  which  have  a  real  edge).  Figure  14(b) 
shows  a  comparison  of  the  standard  MM  and  the  sweep  2  currents  for  the  strip.  Note  that 
the  sweep  2  currents  have  picked  up  the  higher  frequency  oscillations  caused  by  the  first 
reflection  of  the  currents  from  the  right  or  trailing  edge  of  the  strip. 

Currents  for  all  sweeps  k  >  2  are  done  in  precisely  the  same  manner,  except  that  odd 
numbered  sweeps  proceed  from  left  to  right  (p  =  1,2,...,  P),  while  even  numbered  sweeps 
proceed  from  from  right  to  left  ( p  =  P,P  -  1,...,1).  Figure  14(c)  shows  that  the  sweep 
k  =  10  TDMM  currents  are  virtually  identical  to  the  standard  MM  currents. 

For  the  above  strip,  Figure  16  shows  the  RMS  error  in  the  TDMM  currents  versus  the 
sweep  k.  Consider  the  curve  for  Dr  =  2.4A  =  the  width  of  the  R-card  for  the  data  of  Figure 
14.  For  sweep  1  the  RMS  error  (as  compared  to  the  standard  MM  currents)  is  about  0.1, 
and  the  error  drops  essentially  monotonically  to  about  0.001  at  sweep  k  —  16.  For  sweeps 
k  >  16,  the  error  essentially  remains  at  0.001.  Although  this  is  a  very  low  error,  from  a 
theoretical  standpoint,  one  would  want  to  know  why  the  error  does  not  go  to  zero.  To 
answer  this  question,  Figure  16  also  shows  the  RMS  error  for  R-cards  of  width  Dr  =  1.2 A 
and  4.8A.  Note  that  as  Dr  increases,  the  final  RMS  error  decreases.  Thus,  the  use  of  the 
R-cards  on  sweep  1  to  electromagnetically  remove  the  right  edges  of  the  strips  involves  an 
approximation  and  introduces  an  error  which  is  not  removed  by  the  subsequent  sweeps.  As 
the  width  of  the  R-card  increases,  the  approximation  error  decreases,  results  in  a  lower  limit 
on  the  RMS  error. 

d.  Future  Work 

The  initial  results  for  the  TDMM  method  are  encouraging.  To  more  fully  understand  the 
capabilities  and  limitations  of  the  method,  we  intend  to  study  the  following: 
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Figure  15:  Computation  of  the  sweep  2  currents  for  a  PEC  strip. 


Figure  16:  The  RMS  error  for  the  strip  current  versus  the  sweep  k,  and  for  various  values  of 
the  R-card  width. 
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1.  Initial  results  show  the  TDMM  method  to  be  about  10  times  faster  than  the  standard 
MM  for  a  strip  with  N  =  1000  unknowns.  However,  we  need  to  make  a  more  systematic 
study  of  the  TDMM  CPU  time,  and  compare  it  with  the  standard  MM  and  other  “fast” 
methods. 

2.  Apply  the  method  to  a  scatterer  such  as  a  circular  cylinder  which  has  a  shadow  bound¬ 
ary. 

3.  Apply  the  method  to  a  non  smooth  body  with  a  point  of  diffraction  in  the  interior  or 
the  body. 

4.  Show  that,  as  opposed  to  most  other  iterative  methods,  the  currents  for  a  second 
incident  plane  wave  can  be  obtained  with  little  additional  effort. 

5.  Show  that  the  fast  multipole  method  (FMM)  can  be  combined  with  the  TDMM  to 
dramatically  reduce  both  the  CPU  time  and  storage. 

6.  Investigate  the  possibility  of  combining  the  TDMM  with  an  asymptotic  technique  such 
as  the  geometrical  theory  of  diffraction  (GTD). 

3.  Junction  Problems 

a.  Wire/Plate  Junction  Near  an  Edge 

One  strategy  for  modeling  a  geometrically  complex  body  is  to  represent  it  as  an  intercon¬ 
nection  of  simpler  shapes.  For  example,  in  “Electromagnetic  Surface  Patch  Code:  Version 
IV”  (ESP4)  [10]  an  arbitrary  conducting  body  is  modeled  as  an  interconnection  of  thin  wires 
and  polygonal  plates.  It  is  not  enough  to  be  able  to  model  wires  and  plates,  but  one  must 
also  be  able  to  treat  plate/plate  and  wire/plate  junctions.  For  example,  Figure  17  shows 
an  ESP4  model  of  a  monopole  antenna  on  the  underside  of  a  helicopter.  The  helicopter  is 
modeled  as  an  interconnection  of  polygonal  plates,  and  a  special  attachment  mode  is  used  to 
properly  treat  the  wire/plate  junction  [11], [12].  The  attachment  mode  is  valid  only  when  the 
wire/plate  junction  is  greater  than  0.1A  from  the  edge  of  the  attached  plate.  In  constructing 
a  polygonal  plate  model  of  an  arbitrary  body,  such  as  the  helicopter  of  Figure  17,  the  plates 
will  typically  be  of  size  0.1A  to  0.2A.  Thus,  on  a  polygonal  plate  model  of  an  arbitrary  body, 
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Figure  17:  A  monopole  on  a  helicopter  is  modeled  as  an  interconnection  of  wires  and  plates. 

the  wire/plate  junction  will  almost  always  be  less  than  0.1  A  from  an  edge,  and  the  simple 
circular  disk  attachment  mode  is  practically  useless. 

Although  we  worked  on  the  problem  of  a  wire/plate  junction  near  an  edge  about  15  years 
ago,  the  solutions  developed  were  too  complicated  to  be  integrated  into  a  general  purpose 
user  oriented  computer  code  [13], [14].  Also,  they  applied  only  when  the  attachment  point 
was  close  to  a  single  edge.  Thus,  over  the  past  year  we  developed  simper  attachment  modes 
which  can  be  integrated  into  the  ESP4  code.  Further,  they  apply  when  the  attachment  point 
is  close  to  a  few  edges,  as  would  be  the  case  for  a  wire  attached  to  a  polygonal  plate  model 
of  an  arbitrary  body. 
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b.  Strip  on  a  Dielectric  Substrate 

A  new  Version  V  of  the  ESP  code  is  presently  under  development  (not  under  JSEP  support) 
which  will  be  able  to  treat  dielectrics.  In  principle,  ESP5  should  be  able  to  treat  printed 
circuit  antennas,  with  the  substrate  modeled  as  a  thin  dielectric  block,  and  the  printed 
antennas  modeled  by  plates.  However,  the  present  code  is  not  applicable  to  printed  circuit 
antennas  since  the  usual  dielectric  modes  are  not  capable  of  representing  the  large  fields 
at  the  junction  of  the  metallic  strip  and  the  dielectric  substrate.  To  accurately  treat  this 
situation,  we  must  develop  a  strip/dielectric  interface  mode.  This  work  will  begin  in  the 
next  reporting  period. 
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1.  Introduction 

Over  this  past  year,  there  have  been  two  significant  developments  under  this  work  unit,  which 
we  believe  will  have  great  impact  on  the  development  of  finite  methods  for  electromagnetic 
modeling.  The  first  is  the  development  and  analysis  of  a  perfectly  matched  anisotropic 
absorber  for  the  finite  element  method  (FEM).  It  is  perfectly  matched  in  the  sense  that 
a  plane  wave  produces  no  reflections  at  the  air/absorber  interface  for  all  frequencies  and 
incidence  angles.  These  type  of  absorbers  are  used  to  truncate  the  computation  domain 
in  differential-based  numerical  methods  and  are  commonly  called  perfectly  matched  layer 
(PML). 

The  first  PML  has  been  developed  by  Berenger  [1]  for  the  finite  difference  time  domain 
(FDTD)  method.  Chew  and  Weedon  [2]  provide  a  systematic  analysis  of  the  PML  in  terms 
of  the  concept  of  “coordinate  stretching”.  They  demonstrate  that  Berenger ’s  modifications 
to  Maxwell’s  equations  can  be  derived  from  a  more  generalized  form  of  Maxwell’s  equations. 
Recently,  it  has  also  been  discovered  by  Sacks  et.  al.  [3]  that  the  reflectionless  properties  of 
a  material  can  be  achieved  if  the  material  is  assumed  to  be  anisotropic.  Unlike  Berenger’s 
approach,  this  one  does  not  require  a  modification  of  Maxwell’s  equations,  making  it  easier 
to  analyze  in  the  general  framework  of  electromagnetics.  In  this  report,  we  plan  to  review 
the  analysis  for  the  anisotropic  PML  and  demonstrate  its  application  to  FEM.  A  comparison 
of  the  two  PML  formulations  as  well  as  a  more  detailed  discussion  of  the  anisotropic  PML’s 
implementation  into  the  FEM  formulation  is  provided  in  [4]. 
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Another  important  step  is  the  adaptation  of  the  anisotropic  PML  to  the  time  domain. 
In  order  for  any  method  to  be  implemented  in  the  time  domain,  the  method  must  be  dy¬ 
namically  stable  and  formulated  with  a  numerically  stable  scheme  [5].  A  dynamically  stable 
system  is  a  system  that  possesses  the  property  of  remaining  finitely  bounded  over  all  time 
when  excited  by  a  finitely  bounded  source.  Numeric  stability  is  concerned  with  the  behavior 
of  how  a  system  responds  to  a  particular  numerical  implementation  of  the  system.  Very 
few  studies  have  been  performed  on  numerical  stability  in  the  electromagnetics  community. 
Among  these  are  Taflove  [6]  for  the  FDTD  method  and  Lee  and  Sacks  [7]  for  the  Whitney 
elements  time  domain  method. 

Currently,  we  are  not  able  to  come  up  with  a  formulation  of  the  anisotropic  PML  that  is 
dynamically  stable.  To  provide  some  insight  into  our  problem,  we  have  begun  to  study  the 
behavior  of  Berenger’s  equations  in  the  time  domain.  Although  these  equations  have  gained 
widespread  use  within  the  FDTD  community,  no  mathematical  analysis  has  been  done  to 
determine  the  stability  associated  with  the  equations.  We  have  done  such  an  analysis,  and 
it  is  presented  in  [8].  A  review  of  the  numerical  stability  analysis  is  provided  here  for  the 
reader. 

The  last  topic  considers  the  development  of  finite  elements  which  are  specifically  tai¬ 
lored  for  electromagnetics.  Both  nodal  [9]  and  edge  elements  [10],  [11]  produce  solutions 
which  have  significant  numerical  dispersion  error  [12]  because  they  are  based  on  polynomial 
approximations  for  the  field  behavior.  This  means  that  to  maintain  a  given  accuracy,  the 
mesh  density  relative  to  a  wavelength  must  increase  as  the  electrical  size  of  the  computation 
domain  increases.  This  increase  may  significantly  limit  the  usefulness  of  FEM  for  electrically 
large  three-dimensional  geometries.  In  this  report,  we  consider  a  recently  introduced  element 
which  we  call  the  solenoidal  edge  element  [13].  As  its  name  implies,  it  is  edge  based  and 
divergenceless  within  the  element,  and  when  the  order  of  the  element  is  high  enough,  the 
numerical  dispersion  is  negligible. 

At  present,  the  solenoidal  elements  have  been  developed  for  the  two-dimensional  (2-D) 
case  in  which  we  formulate  the  problem  as  a  vector  one  rather  than  a  scalar  one.  We  begin 
with  the  mathematical  analysis  of  the  solenoidal  elements.  The  incorporation  of  higher  order 
absorbing  boundary  conditions  is  especially  easy  with  this  formulation.  We  then  analytically 
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Figure  18:  Plane  wave  incident  on  an  interface  separating  two  halfspaces. 

as  well  as  numerically  show  that  these  elements  produce  very  little  numerical  dispersion. 
Finally,  we  demonstrate  the  effectiveness  of  these  elements  for  scattering  problems. 

2.  Anisotropic  PML  for  Frequency  Domain  FEM 

a.  Review  of  anisotropic  PML 

Referring  to  Figure  18,  the  time-harmonic  form  of  Maxwell’s  equations  in  Region  2  can  be 
written  as 

V  •  [t]E  —  0 

v-[p]ff  =  0 

V  x  E  =  —juj[ji}H 

VxH  =  ju[e\E  (11) 
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where  [/Z]  and  [e]  are  the  effective  permeabihty  and  permittivity  of  Region  2,  respectively. 
In  this  paper,  we  concentrate  on  materials  with  [/Z]  and  [e]  diagonal  in  the  same  coordinate 
system. 
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Furthermore,  we  select  [e]  and  [fi]  such  that 

M  =  M 

e  n 

to  match  Region  2  to  Region  1.  Thus,  we  can  write, 

fa  0  0 
[/Z]  =  /t[A]  =  fi  0  b  0 

V  0  0  c 

fa  0  0 
[e]  =  e[A]  =  e  I  0  b  0 

V  0  0  c 

where  a,b,c  are  in  general  complex.  Consequently,  Eq.  (11)  reduces  to 

V-e[A]i  =  0 
V  ■  fi[A]H  =  0 

Vxi  =  —jufi[A]H 
V  x  H  =  ju)e[A]E 

The  two  curl  equations  in  (14)  can  be  combined  to  form 


(12) 


(13) 


(14) 


[A]  1 V  x  [A]  'Vxi  -  k2E  =  0  (15) 

As  in  Berenger’s  PML,  we  must  specify  a  further  condition  to  the  PML  material  to  make 
it  reflectionless  to  a  plane  wave.  For  the  coordinate  system  defined  in  Figure  18,  we  must 
specify  a  =  b  =  1/c.  For  a  propagating  wave  to  decay  exponentially,  we  must  specify  that 
the  Im(a)  <  0.  Also,  the  Re(o)  can  be  specified  for  the  desired  absorption  of  evanescent 
waves. 
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Figure  19:  Geometry  of  the  PML  region  surrounding  the  scatterer,  regular  view  and  cross 
sectional  view. 

b.  Modeling  of  Edge  and  Corner  Regions  of  the  PML 

The  theory  for  the  anisotropic  PML  is  based  on  the  assumption  that  a  plane  wave  is  prop¬ 
agating  through  a  planar  media  interface  of  infinite  extent.  However,  the  main  interest  for 
the  PML  in  our  case  is  to  use  it  to  absorb  the  scattered  field  from  an  object  in  free  space. 
Thus,  the  PML  material  must  totally  surround  the  scatterer,  and  the  outer  surface  of  the 
PML  must  also  be  terminated  in  some  way.  For  this  work,  we  terminate  the  PML  with  a 
perfect  electric  conductor.  The  PML  material  is  placed  in  the  shape  of  a  box  (Figure  19)  in 
order  to  best  approximate  the  reflectionless  properties  of  the  PML. 

The  choice  for  the  material  properties  of  the  side  regions  of  the  box  (A^lde,  i  =  x,y,z)  is 
straightforward  and  is  shown  in  the  figure,  where  sx,sy,sz  are  arbitrary  complex  numbers. 
However,  the  method  for  determining  material  properties  at  the  edge  and  corner  regions  of 
the  box  is  not  so  clear.  From  electromagnetic  theory,  one  expects  that  any  choice  would 
cause  some  diffracted  field  from  the  edges  and  corners.  One  approximate  approach  for  the 
edge  region  is  to  choose  the  edge  properties  such  that  they  are  perfectly  matched  to  the 
adjacent  side  regions  when  the  edge/side  interface  is  of  infinite  extent.  We  can  use  a  similar 
approach  for  the  corner  region  by  matching  the  corner  properties  to  the  adjacent  edges.  If 
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Table  1:  Parameters  for  the  four  test  cases  used  in  rectangular  plate  example. 


di(  A) 

<*2(A) 

a 

P 

meshsize  (A) 

no.  of  unknowns 

Case  1 

0.2 

0.2 

1.5 

1.5 

0.05 

212976 

Case  2 

0.1 

0.2 

1.5 

1.5 

0.05 

129396 

Case  3 

0.1 

0.2 

2.0 

2.0 

0.05 

129396 

Case  4 

0.2 

0.2 

4.0 

4.0 

0.05 

212976 

we  go  through  the  corresponding  analysis  of  a  plane  wave  propagating  through  an  infinite 
media  interface  for  the  edges  and  corners,  we  arrive  at  the  interesting  relationship  that  A edge 
is  equal  to  the  matrix  multiplication  of  A S1(fe  of  the  two  adjacent  sides  and  Acorner  is  equal  to 
the  matrix  multiplication  of  A side  of  the  three  adjacent  sides.  As  an  example,  let  us  consider 
the  edge  in  Figure  19  with  material  properties  Kvedge.  The  matching  condition  is, 


t  SyS2 
0 
0 


0  0  \ 

Sz/Sy  0 

0  Sy/Sz 


c.  Numerical  results  for  anisotropic  PML 


(16) 


In  order  to  test  the  capabilities  of  the  anisotropic  PML,  we  consider  its  use  in  the  edge- 
based  finite  element  method  [10].  The  implementation  of  the  anisotropic  media  into  the 
finite  element  formulation  is  simple  and  straightforward  [14]  and  therefore  is  not  discussed 
here.  A  bi-conjugate  gradient  solver  (BiCG)  with  diagonal  preconditioning  is  used  to  solve 
the  matrix  equation.  The  material  properties  for  all  the  numerical  results  are  given  by 
sx  =  sy  =  sz  =  a  —  j(3,  where  sx,  sy,  sz  are  as  shown  in  Figure  19,  and  a  and  (3  are  some 
real  numbers. 

For  the  first  example,  we  consider  the  scattering  of  plane  waves  by  a  larger  rectangular 
plate  in  Figure  20.  The  size  of  the  plate  is  2A  x  1A,  and  the  mesh  is  uniform.  Our  purpose  is 
to  try  different  values  of  <fj,  d2>  a,  and  /3,  and  see  how  they  affect  the  patterns  for  this  plate. 

Figure  21  shows  the  ECS  patterns  for  the  four  different  cases  listed  in  Table  1  as  compared 
to  the  method  of  moments  (MOM)  results.  As  seen  in  Table  1,  Cases  1  and  2  have  the  same 
material  property  but  different  Cases  2  and  3  have  the  same  mesh  but  different  material 
properties;  Cases  1  and  4  also  have  the  same  mesh  but  different  material  properties.  As  seen 
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di 


di 


in  Figure  21,  Case  1  provides  better  absorption  than  Case  2;  however,  when  larger  material 
property  than  those  for  Case  2  is  used,  such  as  in  Case  3,  the  absorption  is  compensated, 
and  very  similar  pattern  as  that  of  Case  1  is  obtained.  Finally,  Case  4  shows  the  pattern 
calculated  by  using  a  very  large  material  property  (a  =  /?  =  4).  It  is  clear  that  the  pattern 
for  Case  4  is  much  worse  than  that  for  Case  1.  This  is  caused  by  the  increases  in  the 
discretization  errors  and  possibly  the  diffractions  by  the  material  interfaces. 

For  the  second  example,  we  consider  the  scattering  of  a  plane  wave  by  a  wedge- cylinder 
plate.  The  problem  geometry  is  shown  in  Figure  22.  The  plate  is  composed  of  a  half  circle 
of  radius  1A  and  an  equilateral  triangle  of  side  length  2A.  The  PML  material  is  placed  as 
close  as  0.1A  to  the  plate  in  the  *  and  y  directions,  and  0.2A  in  the  z  direction. 

It  should  be  mentioned  that  a  non-uniform  mesh  is  used  for  this  example.  The  mesh  is 
generated  in  such  a  way  that  it  is  denser  inside  the  PML  material,  in  the  directions  normal  to 
the  free-space/PML  interfaces,  where  the  fields  would  change  the  most  rapidly.  It  is  hoped 
that  by  so  doing  the  number  of  unknowns  can  be  reduced  without  too  much  increase  in  the 
discretization  errors.  The  resulting  number  of  unknowns  is  48242,  in  a  problem  domain  of 
volume  6.968A3,  as  compared  to  129396  out  of  a  volume  of  2.496A3  for  Case  2  of  the  previous 
example.  So  the  average  mesh  size  is  much  larger  than  the  previous  example. 

Figures  23  and  24  show  the  bistatic  RCS  patterns  ( x-z  plane;  4>  polarization)  for  two 
incident  waves  6inc  =  <f>tnc  =  0°  and  6'nc  =  30°,  <f>mc  =  0°,  respectively,  both  having  y- 
polarized  electric  field.  The  material  property  a  =  (3  =  1.5  is  used  for  both  cases.  It  is  seen 
that  excellent  agreement  with  the  MoM  patterns  are  obtained  for  both  cases. 

3.  Numerical  Stability  Analysis  for  Berenger’s  PML 

Assuming  for  the  moment  that  Berenger’s  equations  are  dynamically  stable,  the  numeric 
stability  of  a  set  of  time  domain  equations  is  examined.  Berenger’s  time  domain  equations 
are  given  as 
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Figure  23:  Bistatic  patterns  for  the  wedge-cylinder  plate.  Patterns  are  in  the  x-z  plane  for  the 
4>  polarization.  The  incident  plane  wave  has  y-polarized  electric  field  with  0tnc  =  <f>inc  =  0°. 
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Figure  24:  Bistatic  patterns  for  the  wedge-cylinder  plate.  Patterns  are  in  the  x-z  plane  for 

the  (f>  polarization.  The  incident  plane  wave  has  ^-polarized  electric  field  with  (j)inc  =  0°  and 
Qinc  _  30°. 
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(17) 


+  £* Hsz  =  —§^z  x  E  ea-§?  +  crzEaz  =  |u5 

Both  the  electric  and  magnetic  total  fields  ( E ,  if)  are  formed  as  the  sum  of  the  respective 
coordinate  stretched  vector  fields.  These  are  expressed  as 

E  =  Eax  +  Eav  +  C  H  =  Hax  +  +  i4  (18) 


while  the  cartesian  total  field  components  are  found  as 

Ex  =  (Eay  +  E„ )  • x  Hx=  (lf3y  -f  Haz)  •  x 

Ey  =  Iesx  +  Eaz)  •  y  Hy  =  \Hax  +  Haz\  •  y  (19) 

Ez  =  yEax  +  Eaf/j  •  z  Hz  =  \Hax  +  HayJ  •  z 

For  convenience,  the  cartesian  components  of  the  coordinate  stretched  vectors  are  ex¬ 
pressed  as 


Tsa0=ls0-CX 

where  T  may  be  either  the  E  or  H  field  and  a ,  /3,  and  (  are  each  uniquely  assigned  to  one 
of  the  cartesian  directions  x,  y,  or  z. 

We  form  twelve  FDTD  equations  using  average  discrete  time  values  where  appropriate 
and  applying  the  central  difference  approximation  to  (17).  By  using  a  leapfrog  scheme  in 
time  and  the  standard  Yee  cells  spatial  distribution,  we  obtain  a  general  formula  for  the  E 
and  H  field  components  as  follows, 


Es^(hhk)  =  f3Esaa(hhk)  +  s9n{oc  x  /3)JJ^ 


+  Sx0,j  +  Sy0,  k  +  6Z0)  -  k ) 


HnaJ{i,i,  *)  =  ffiKJfrj, k )  +  39n(a  x 

[Ec  (i  —  bxfiij  ^y/3?  k  &zfi)  “1"  (*'»  Jj  ^)] 


where  fa  and  ja  are  defined  as 


7Q  =  2e  +  cra  At  fa  = 


2e— <Tg  At 


and  Sa0  is  defined  as 


63 


&a0  — 


if  a  and  /3  are  the  same 
cartesian  components 
otherwise 


(23) 


Note  that  as  the  range  of  aa  is  varied  from  0  to  oo,  fa  spans  from  1  to  -1.  Next,  the  E 
field  components  of  (21)  are  expanded  so  that  they  are  only  a  function  of  the  n  and  n—  |  time 
steps.  In  this  way  both  field  components  have  the  same  forms.  Combining  Equations  (19) 
through  (21)  yields 
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As  an  example,  the  Esyx  component  of  these  modified  equations  is  worked  out  explicitly. 
For  this  component  the  E  field  equations  are  used  while  a  =  y,  (3  =  z,  and  (  is  selected  as 
the  z  cartesian  component. 


i,k) 


If »  (#.”/(•  +  1  ,j,k)  -  ))  + 

t,  + 1  ,j,k)  -  b&v,  j,k)) }  + 

K  (»  +  1J-  1,  k)  +  (i,i,  k) 


e4A  e 

fiAyAx'fy'/x 


-E:x(i,j-i,k)-E:x(i+ij,k)] 

e  (  2A t  \2r 

--  j  ^ + u,*) 

1  ,},k)] 


(25) 
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Following  the  Von  Neumann  analysis  the  electric  and  magnetic  fields  are  separated  into 
time  varying  modal  series, 

En(i,j,  k)  =  E  Y,  E  en(t)e5(miAx+njAy+pkAz)  (26) 

m  n  p 

Hn+*(i,j,k )  =  EE  ^hB+?(t)ei(miAx+njAy+pkAz)  (27) 

m  n  p 

where  j  =  y/—l  and  e(t)  and  h(t)  are  complex  functions  of  time.  Substituting  these  series 
into  (25),  we  obtain  a  matrix  equation  in  terms  of  the  matrix  [G], 


Un+1  =  [G]Un 


where 
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[G22]  — 


(34) 


'  o  0  0  0  O' 

0  /*  0  0  0  0 

0  0  /„  0  0  0 

0  0  0  /„  0  0 

0  0  0  0  fz  0 

.0  0  0  0  0  fz  . 

where 

Cx  =  -e-jiAx  Cy  =  Cz  =  -e~ikAz  (35) 

/*  "  V 

k-  =  ~  V  k*  =  2^’  -  1)  *•  =  -  1)  (36) 

Stability  is  now  reduced  to  finding  the  eigenvalues  of  [G]  and  showing  the  conditions  under 
which  the  modulus  of  these  eigenvalues  are  bounded  by  one.  The  eigenvalues  of  [G]  are  found 
from  the  solution  of  the  characteristics  equation.  This  characteristic  equation  decomposes 
into  a  repeated  sixth  order  polynomial  equation.  This  implies  that  all  eigenvalues  A  are 
repeated.  Finding  the  six  unique  eigenvalues  reduces  to  finding  the  zeroes  of  the  following 
sixth  order  equation  (5(A), 

<5(A)  =  A  +  a5A5  +  CI4A4  +  a3A3  +  a^A2  +  01A  +  ao  (37) 

where 

a5  =  -(X  +  Y  +  Z  +  2[fx  +  fy  -f  /,])  (38) 

=  2  (X[fy  +  fz]  +  Y[fx  +  fz]  +  Z[fx  +  fy])  + 

4(y tf y  +  fxfz  +  fyfz)  +  fx+  fy  +  fl  (39) 

«3  =  -m2y+fl)  +  Y(fx+fl)  +  Z(fx+fy)  + 

2  (f2Afy  +  /*]  +  Pyifx  +  fz]  +  f2z[fx  +  fy])  + 

W  +  Yfxfz  +  Zfxfy)  +  8  (fjyfz))  (40) 

02  =  2{X[fyfZ  +  flfy]  +  Y[flfZ  +  flfX]  +  Z{flfy  +  f2yfX])  + 
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(41) 


4  (fxfyfz  +  fxfyfz  +  fxfyf)  +  ffy  +  fxfz  +  fyfz 

*1  =  -{Xflf2z+YflJl  +  Z}lfl  + 

2  {fxflfl  +  flfyfl  +  fxfyfz))  (42) 

OO  =  fxfyfz  (43) 

and 
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P  ' 

[a z7z)  2 

Note  that  because  X,  Y ,  and  Z  are  purely  real,  the  coefficients  of  (37)  are  also  purely 
real. 

Next  a  simple  case  is  studied  and  a  generalized  stability  criteria  follows.  In  this  example 
case  the  wave  is  absorbed  in  only  one  direction  (say  the  z-direction)  which  implies  that 
(< crx  =  <Ty  =  0)  and  (erz  ^  0).  Further,  let  the  spatial  resolution  be  the  same  in  all  directions 
(A®  =  Ay  =  A z).  This  implies  that  (fx  =  fy  =  1)  while  ( fz  1). 

For  stability,  the  values  of  X,  Y ,  and  Z  are  chosen  so  as  to  produce  the  eigenvalues  with 
the  largest  moduli.  The  stability  is  therefore  satisfied  for  all  other  choices  of  indexes  i,  j, 
and  k. 

Putting  this  information  into  (37),  it  is  obvious  that  (A  —  l)2  and  P( A)  both  divide  £?(A), 
where  P( A)  is  the  fourth  order  polynomial 

P( A)  =  04A4  +  a3A3  -j-  <i2^2  -h  aiA  -(-  do  (47) 

with  coefficients 


a4  =  1 

(48) 

a3=  2(2p2(2+(^f±) 

2)  -(/.  +  !)) 

(49) 

a2=  _(8p2(2/,  +  (^±l)2) 

-  (fz  +  4 fz  +  1)) 

(50) 
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where 


ai  = 


(51) 

(52) 


2(2^(2/,2  +  (If l)2)  -  Mf,  +  1)) 

“o  =  /,2 


Ai 

fieAx 


(53) 


Using  the  Jury  test  [15],  the  modulus  of  the  roots  are  forced  less  than  one  if  the  following 
5  conditions  are  met:  P(l)  >  0,  P(-l)  >  0,  \a0\  <  a4?  |50|  >  |63|,  and  |c0|  >  |c2|  where  bk 
and  ck  are  defined  as 


do 

q<4 

b0  63- 

&3 


The  above  five  conditions  are  all  satisfied  when 


(54) 


P2  <  1  (55) 

Note  that  the  numerical  stability  of  Berenger’s  formulation  is  not  dependent  on  the 
choice  of  az  >  0.  For  the  general  case,  the  values  of  ax  >  0,  Cy  >  0,  and  az  >  0  can  be 
chosen  without  effecting  the  numerical  stability.  Even  in  the  limit  as  all  three  of  these  values 
approach  infinity,  Equation  (37)  reduces  to 

Q(X)  =  (A  +  l)6  (56) 

and  the  roots  remain  bounded  by  1. 

4.  Solenodial  Edge  Elements 

a.  Analysis 

The  solenoidal  edge  elements  have  been  developed  specifically  to  minimize  error  due  to 
numerical  dispersion.  These  elements  must  satisfy  two  criteria  other  than  those  imposed  by 
good  mesh  generation  practices.  Namely,  all  elements  must  be  triangular,  and  the  material 
must  be  homogeneous  within  each  element.  A  short  review  of  the  solenoidal  edge  elements 
is  presented  in  this  section. 
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Figure  25:  Cylindrical  coordinate  system  for  definition  of  the  basis  functions  associated  with 
the  edge  number  1. 

The  unknown  field  is  expanded  in  terms  of  basis  functions,  separately  within  each  element. 
These  are  vector  basis  functions,  and  as  its  name  implies,  they  are  edge  based  and  divergence 
free  within  the  element.  We  require  that  each  basis  function  within  an  element  be  associated 
with  only  one  of  the  edges.  This  can  be  easily  achieved  by  forcing  tangential  components  of 
the  vector  basis  to  zero  along  all  of  the  remaining  edges. 

To  be  more  specific,  let  us  analyze  the  basis  function  associated  with  edge  number  1  of 
an  isolated  triangular  element  as  shown  in  Figure  25.  We  us  assume  that  the  polarization 
is  T Ez  although  the  derivation  for  the  T Mz  case  follows  similarly.  According  to  the  above 
discussion,  the  vector  basis  must  have  zero  tangential  components  along  edges  2  and  3,  and 
since  the  transversal  field  that  we  are  presently  analyzing  is  an  electric  field,  it  is  immediately 
seen  that  the  same  result  is  obtained  if  the  edges  2  and  3  are  made  of  perfectly  conducting 
material. 
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We  can  now  visualize  these  perfectly  conducting  edges  as  forming  the  tip  of  an  infinitely 
long  wedge  with  perfectly  conducting  walls,  as  in  Figure  26.  An  analysis  of  the  fields  satisfy¬ 
ing  the  wave  equation  inside  such  an  infinite  wedge  shows  that  the  wedge  supports  infinitely 
many  standing  wave  modes.  It  is  important  to  note  that  the  order  of  the  mode  is  determined 
by  an  integer  n  and  that  the  higher  order  modes  exhibit  faster  spatial  variation  than  the 
lower  order  modes.  These  characteristics  of  the  modes  are  identical  to  those  that  we  require 
of  an  edge  based  solenoidal  vector  basis.  The  solenoidal  edge  basis  for  edge  number  1  in 
Figure  25  is  defined  as: 

/  i\  Jisnikp)  .  , 

WPn[P^(P)  =  Vn - Yp - Sin\Vn9) 

™4>n{p,<t>)  =  j'vn{kp)cOs(l/n<f>)  (57) 

K{p,<!>)  =  ~curl(wn)  ■  z  =  -J^kpyos^cf)) 

7T 

where:  vn  =  n— ,  n  =  0,1,2,... 

0m 


and  k  and  Z  are  the  wavenumber  and  the  characteristic  impedance  for  the  homogeneous 
material  inside  the  element.  The  hn  basis  is  used  for  the  Hz  field  expansion  and  relates 
directly  to  the  curl  of  the  wn  basis  for  the  transversal  E  field.  A  different  cylindrical  coor¬ 
dinate  system  must  be  used  for  each  edge;  the  center  of  the  coordinate  system  is  always  at 
the  node  opposing  the  edge  in  question. 

Let  us  consider  some  notation  for  the  next  equation.  We  use  index  e  to  denote  the  element 
number,  index  i  for  the  edge  number,  and  index  n  for  the  order  of  the  basis  function.  We 
assume  that  the  highest  order  of  the  basis  for  the  ith  edge  is  given  by  jV,-.  The  electric  and 
magnetic  field  of  a  TEZ  polarized  wave  can  now  be  approximated  by: 


Ni-l 


=  EEE«'W 


i  n— 0 

=  fEEEWB  (58) 

where  the  first  summation  is  over  all  elements,  the  second  summation  is  over  all  three  edges 
belonging  to  the  element,  and  the  third  summation  is  over  TV,  orders  of  the  basis  function. 
Basis  function  w ®’*  is  the  nth  order  basis  associated  with  the  eth  element  and  edge  i  which 
belongs  to  that  element.  The  coefficients  are  the  amplitudes  of  the  basis  functions  in 
the  expansion  for  the  field,  and  they  are  therefore  the  unknowns  in  our  model. 
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Figure  26:  Infinite  corner  with  perfectly  conducting  walls  and  filled  with  a  homogeneous 
material  of  characteristic  impedance  Z  and  wavenumber  k. 
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We  see  from  (58)  that  the  wave  equation  is  exactly  satisfied  within  each  element  even 
before  we  proceed  with  the  computation  of  the  unknown  coefficients  This  is  the  reason 
that  very  little  numerical  dispersion  is  present  in  the  numerical  solution.  To  find  the  unknown 
coefficients,  we  apply  continuity  of  the  tangential  electric  and  magnetic  fields  across  the 
interelement  boundaries.  Let  us  focus  our  attention  to  the  problem  of  field  matching  at  the 
interelement  boundaries,  which  is  indeed  the  key  step  in  the  implementation  of  the  solenoidal 
edge  elements.  First,  we  should  realize  that  each  edge  in  the  interior  of  the  mesh  belongs 
to  two  elements,  both  of  which  have  separate  basis  functions  associated  with  this  edge.  We 
require  that  the  number  of  these  bases  on  both  sides  of  the  edge  be  equal.  Let  this  number 
of  basis  functions  associated  with  the  edge  i  be  equal  to  JV,-.  For  each  interior  edge  i  there 
are  two  elements  e\  and  e 2  which  share  the  edge.  The  number  of  unknowns  associated  with 
the  edge  is  simply  twice  the  number  of  bases  on  each  side,  that  is  2 IV,-.  Thus,  the  matching 
condition  for  edge  i  must  be  expressed  in  terms  of  21V,  equations.  We  form  IV,  equations  by 
matching  the  tangential  electric  field  and  the  other  IV,  equations  by  matching  the  magnetic 
field  along  the  ilh  edge. 

The  field  matching  is  performed  in  an  integral  sense.  The  matching  equations  for  the 
electric  field  defined  by  the  edge  common  to  elements  e1  and  e2  can  be  put  down  in  the 
following  form: 

J0  n=0 

=  All  <”<”(£.■)] 

J0  n= 0 

for:  p  =  -  1  (59) 

where  t,  is  the  unit  vector  of  the  ith  edge.  The  weighting  function  is  given  by 

*»(&)  =  «w(»t6),  0  <  £,<  1,  n  =  0, ...,  Ni  —  1  (60) 

where  £,  is  the  normalized  linear  coordinate  along  the  ith  edge.  Matching  of  the  magnetic 
field  can  be  done  in  the  same  manner,  with  an  important  difference  that  the  tangential 
magnetic  field  along  the  edge  z,  depends  not  only  on  basis  functions  associated  with  this 
edge,  but  also  on  all  other  basis  functions  defined  within  the  elements  e a  and  e2.  It  should 
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be  noted  that  ABC’s  are  easily  implemented  into  the  solenoidal  element  formulation  because 
the  derivatives  of  the  basis  functions  can  be  evaluated  analytically. 

b.  Numerical  results  for  solenoidal  elements 

For  the  scattering  case,  our  main  interest  is  in  the  interaction  of  the  ABC  with  the  solenoidal 
elements.  However,  before  we  do  that,  let  us  consider  the  issue  of  numerical  dispersion.  To 
validate  the  claims  of  low  numerical  dispersion,  we  consider  the  problem  of  plane  wave  ( TMZ 
polarization)  scattering  from  a  homogeneous  dielectric  circular  cylinder  of  radius  R\  —  0.8Ao 
(A0  is  the  free  space  wavelength)  with  erj  =  4.  Since  the  exact  solution  is  known,  we  can 
quantify  the  error  by  defining  the  %  RMS  error  as 

%  RMS  Error  =  j ^  *  X  WO  (61) 

where  E^x  is  electric  field  produced  by  the  exact  eigenfunction  solution  with  the  radiation 
boundary  condition  and  E^BC  is  the  electric  field  calculated  with  a  first  order  BGT  ABC  as 
the  outer  boundary  condition.  The  term  E^BC  represents  the  solution  due  to  either  bilinear 
quadrilateral  nodal  elements,  the  solenoidal  elements,  or  the  analytical  solution  with  the 
ABC  [16].  The  integration  path  Tj  is  chosen  to  be  the  along  the  surface  of  the  cylinder. 

The  %  RMS  error  is  shown  in  Figure  27  as  R2,  which  represents  the  location  of  the  outer 
boundary,  is  varied  from  0.9Ao  to  2Ao  for  the  analytical,  solenoidal  elements  (mesh  spacing 
of  A0/15  with  third  order  elements),  and  nodal  elements  (mesh  spacings  of  approximately 
Ao/20  and  Ao/40).  The  error  in  the  analytical  solution  is  due  solely  to  the  ABC.  The  curve 
of  the  error  in  the  solenoidal  element  follows  very  closely  to  the  error  curve  for  the  analytical 
case,  which  indicates  that  the  solenoidal  elements  produce  very  little  numerical  dispersion. 
The  slight  differences  are  most  likely  caused  by  the  piecewise  linear  approximation  of  the 
circular  shape  of  the  cylinder  which  is  used  by  the  mesh  for  the  solenoidal  elements.  On  the 
other  hand,  the  deviation  of  the  nodal  element  errors  from  the  analytical  case  is  due  almost 
entirely  to  numerical  dispersion. 

To  study  the  accuracy  of  the  solenoidal  elements  with  the  ABC,  we  consider  the  problem 
of  a  TEZ  polarized  plane  wave  incident  on  a  dielectric  circular  cylinder.  The  geometry  is 
shown  in  Figure  28,  and  the  associated  mesh  is  given  in  Figure  29  with  282  elements  and 
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Figure  27:  Plot  of  %  RMS  error  on  the  dielectric  circular  cylinder  boundary  as  function 
of  R2/A0  for  the  analytical,  solenoidal  solution,  and  nodal  (Ao/20  and  Ao/40  mesh  spacing) 
cases,  i?!  =  0.8Ao,  erd  =  4. 
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Figure  28:  Geometry  for  a  radiation  type  test  problem.  Dielectric  cylinder  with  er  =  2  and 
radius  0.4A  is  illuminated  by  a  TE  polarized  plane  wave.  Fields  will  be  compared  along  the 
sampling  line. 


Figure  29:  Mesh  for  a  radiation  test  problem  shown  in  Fig.  26.  An  absorbing  boundary 
condition  is  applied  on  the  square  boundary  of  the  mesh. 
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439  edges.  In  the  following  plots,  we  look  at  the  electric  field  along  the  sampling  line  shown 
in  Figure  28.  For  the  first  case,  we  use  first  order  elements  with  a  first  order  ABC.  The 
magnitude  of  the  x  and  y  components  of  the  electric  field  is  plotted  in  Figure  30  for  the 
solenoidal  element  and  compared  to  the  series  solution.  The  two  components  in  the  FEM 
solution  behave  differently  because  the  mesh  is  oriented  such  that  Ey  is  tangential  to  the 
edges  while  Ex  is  normal  to  the  edges.  The  first  order  basis  function  produces  an  almost 
constant  tangential  field  along  the  edge  and  a  linear  behavior  for  a  field  normal  to  the  edge. 
We  next  consider  the  case  where  first  order  elements  are  used,  but  a  second  order  ABC  is 
placed  at  the  boundary.  In  Figure  31,  we  plot  the  magnitude  of  the  differences  between  the 
numerical  and  series  solution  electric  field  for  first  order  solenoidal  elements  with  either  a 
first  or  second  order  ABC.  Interestingly  enough,  the  accuracy  of  the  solution  is  worse  for  the 
higher  order  ABC.  It  is  likely  that  much  of  the  error  is  due  to  poor  discretization.  Also,  the 
failure  of  the  second  order  ABC  to  outperform  the  first  order  one  may  very  likely  be  due  to 
the  inaccuracies  in  the  first  order  solenoidal  elements.  In  order  to  test  this  hypothesis,  we  use 
third  order  solenoidal  elements  with  a  first  order  ABC  on  the  same  geometry  as  Figure  28 
with  the  same  mesh.  The  electric  field  components  are  plotted  in  Figure  32.  Note  that  the 
solution  is  much  smoother  and  much  more  accurate  than  for  the  first  order  elements.  If  we 
now  use  a  second  order  ABC,  we  obtain  much  more  accurate  results  (Figure  33)  than  for  the 
first  order  ABC.  These  results  demonstrate  that  in  order  to  gain  the  advantages  of  a  higher 
order  ABC,  we  must  choose  the  element  order  to  be  high  enough. 

To  demonstrate  the  accuracy  of  these  elements  for  an  electrically  larger  case,  we  in¬ 
crease  the  electrical  size  of  the  geometry  in  Figure  28  by  a  factor  of  five.  Thus  the  radius 
of  the  dielectric  circular  cylinder  is  2Ao.  Using  the  same  mesh  as  in  Figure  29  (1.6  ele¬ 
ments/wavelength),  we  can  produce  a  very  accurate  solution  (Figure  34)  with  third  order 
elements  and  a  second  order  ABC. 
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Figure  31:  Plot  of  the  magnitude  of  the  error  for  the  transversal  electric  field  for  the  case  in 
Figs.  27  and  28.  Errors  are  compared  for  the  solenoidal  edge  elements  that  use  1st  and  2nd 
order  ABC.  1st  order  elements  are  used  in  both  cases. 
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Figure  32:  Plot  of  the  magnitude  of  the  vertical  and  horizontal  electric  field  components  for 
the  case  in  Figs.  27  and  28.  3rd  order  solenoidal  edge  element  solution  with  1st  order  ABC 
is  compared  to  the  exact  solution. 
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Figure  33:  Plot  of  the  magnitude  of  the  error  for  the  transversal  electric  field  for  the  case  in 
Figs.  27  and  28.  Errors  are  compared  for  the  solenoida!  edge  elements  that  use  l4t  and  2nd 
order  ABC.  3rd  order  elements  are  used  in  both  cases. 


Vertical  position  y 


Figure  34:  Plot  of  the  magnitude  of  the  vertical  and  horizontal  electric  field  components  for 
the  case  in  Figs.  27  and  28,  but  for  5  times  electrically  larger  problem.  3rd  order  solenoidal 
edge  element  solution  with  2nd  order  ABC  is  compared  to  the  exact  solution. 
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1.  Introduction 


EM  phenomena  associated  with  radiating  objects  can  be  classified  according  to  different 
scales  in  terms  of  the  radiated  EM  wavelength.  For  example,  at  one  extreme,  for  the  high 
frequency  regimes  of  interest,  a  complex  radiating  object  such  as  an  aircraft  can  exhibit  both 
electrically  very  large  as  well  as  electrically  small  portions;  at  the  other  extreme,  radiating 
antenna  elements  and  their  feed  components  can,  either  due  to  need  for  miniaturization  or 
due  to  use  of  very  high  frequencies,  become  extremely  small  physically  (as,  for  example, 
in  some  modern  wireless  systems,  or  as  in  millimeter  wave  antenna  array  elements,  etc.) 
requiring  one  to  consider  integrating  a  small  radiating  element  directly  with  its  source  which 
is  a  far  smaller  active  solid  state  device.  Both  types  of  the  above  situations,  depicting 
extremes  of  scales,  are  expected  to  occur  even  more  frequently  in  the  future.  However,  there 
is  no  single  approach  or  technique  which  can  be  used  as  a  prediction  or  analysis/design  tool 
for  dealing  with  such  EM  phenomena  in  a  tractable  fashion.  For  example,  neither  asymptotic 
high  frequency  techniques,  nor  numerical  or  the  low  frequency  techniques  alone  can  predict 
the  radiation  or  scattering  from  an  object  containing  large  and  small  parts;  also  full  wave 
techniques  which  can  analyze  an  antenna  element  cannot  really  easily  deal  with  an  active 
solid  state  device  which  may  be  used  as  a  source  for  the  antenna.  It  thus  becomes  extremely 
important  to  develop  systematic  hybrid  combinations  of  different  (two  or  more)  techniques 
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in  which  each  technique  is  respectively  chosen  so  that  it  is  best  suited  for  efficiently  handling 
a  different  part  of  the  complete  radiating  object,  and  these  techniques  when  systematically 
combined  in  a  given  hybrid  combination  can  then  handle  the  entire  or  complete  radiating 
object  in  a  tractable  fashion.  Clearly,  the  development  of  such  hybrid  tools  would  be  of 
significant  interest,  and  also  of  use  to  present  and  future  EM  technology.  With  this  view  in 
mind,  this  hybrid  section  of  JSEP  has  focused  on  the  development  of  hybrid  prediction  or 
analysis/design  tools. 

During  the  past  period,  substantial  progress  was  made  in  the  development  of  a  hybrid 
method  for  analysis/design  of  an  active  integrated  microstrip  patch  antenna  using  a  combi¬ 
nation  of  circuit  solvers  and  nonlinear  device  analysis  to  deal  with  an  active  bipolar  junction 
transistor  (BJT)  operating  as  an  oscillator,  with  a  full  wave  analysis  for  the  microstrip  patch 
type  antenna  radiating  element  whose  parameters  were  integrated  into  the  design  of  the  ac¬ 
tive  device  operation.  This  active  integrated  antenna  is  an  example  where  it  is  not  possible 
to  design  the  passive  radiating  patch  with  some  idealized  feed  as  is  often  done  in  microstrip 
patch  antenna  full  wave  analysis;  indeed,  it  is  important  to  know  the  actual  active  device 
parameters  as  closely  as  possible  so  that  the  passive  resonant  load  provided  by  the  patch  can 
be  truly  integrated  into  the  active  device  analysis  so  that  the  active  device  can  provide  stable 
oscillations  at  the  desired  frequency  of  operation.  A  careful  large  signal  analysis  of  the  BJT 
as  well  as  an  accurate  analysis  of  the  antenna  is  required  to  conduct  a  modified  Kurokawa’s 
nonlinear  analysis  to  guarantee  stable  oscillations  within  its  active  region.  Note  that  this 
analysis  takes  into  account  several  harmonics  in  addition  to  the  fundamental  signal. 

Also,  in  the  development  of  hybrid  methods,  it  is  necessary  to  consider  the  use  of  partial 
differential  equation  (PDE)  based  numerical  methods  for  solving  Maxwell’s  equations  which 
are  especially  useful  in  dealing  with  regions  of  a  radiating  object  which  are  electrically 
small  but  highly  inhomogeneous.  The  finite  different  time  domain  (FDTD)  has  become 
a  rather  popular  numerical  PDE  technique  despite  some  drawbacks  such  as,  for  example, 
its  difficulties  in  handling  curvilinear  geometries;  on  the  other  hand,  the  computational 
fluid  dynamics  (CFD)  community  appears  to  have  made  significant  strides  in  PDE  based 
numerical  methods  which  has  led  to  the  finite  volume  time  domain  (FVTD)  approach  that 
appears  to  be  potentially  more  advantageous  especially  from  the  point  of  view  of  its  versatility 
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in  dealing  with  curvilinear  boundaries,  and  its  robustness.  However,  a  thorough  comparison 
of  FDTD  with  FVTD  does  not  appear  to  have  been  conducted  making  it  unclear  as  to 
which  method  is  better  suited  for  what  situations.  Consequently,  a  study  of  one  form  of  this 
potentially  useful  FVTD  has  begun  to  characterize  its  efficiency,  stability,  use  in  wide  band 
and  even  non-liner  regimes,  and  its  behavior  with  respect  to  absorbing  boundary  conditions 
for  truncation  of  FVTD  regions.  The  initial  studies  are  focusing  on  antenna  phenomena  in 
the  presence  of  antenna  excitation,  and  in  presence  of  materials  for  optimizing  the  antenna 
performance.  An  elliptic  body  fitting  grid  generator  has  also  been  developed  in  2-D  to  be 
used  in  the  FVTD  discretization.  It  is  hoped  that  this  study  would  provide  more  useful 
information  on  the  FVTD  and  to  also  indicate  areas  of  further  study  to  possibly  improve 
the  efficiency  of  the  FVTD  procedure  while  still  maintaining  its  versatility. 

Work  has  also  continued  on  the  development  of  a  hybrid  method  for  dealing  with  material 
treatment  of  edged  bodies  in  which  the  material  can  control  the  scattering  characteristics 
of  the  edge.  An  area  where  this  hybrid  approach  has  an  impact  is  in  antenna  design.  Most 
antenna  structures  contain  edges  and  the  treatment  of  edge  effects  becomes  crucial  in  many 
applications  to  obtain  the  desired  antenna  performance.  The  use  of  the  special  Green’s  func¬ 
tion  for  the  impedance  wedge  developed  in  the  Diffraction  Studies  unit  has  been  employed 
in  this  hybrid  approach  to  demonstrate  its  usefulness  in  the  design  of  a  microstrip  antenna 
where  edge  effects  produce  detrimental  effects  on  the  antenna  performance,  and  where  these 
undesired  effects  are  nullified  to  obtain  an  excellent  antenna  performance  by  proper  design 
of  edge  treatment  using  a  hybrid  numerical  moment  method  asymptotic  approach.  In  this 
hybrid  approach,  an  integral  equation  formulation  is  employed  which  is  solved  by  a  numer¬ 
ical  moment  method  (MM)  procedure,  but  the  number  of  unknowns  to  be  solved  in  this 
MM  is  drastically  reduced  because  of  the  special  choice  of  the  aforementioned  wedge  Green’s 
function  in  the  kernel  of  the  integral  equation. 

Significant  progress  has  also  been  made  in  the  hybrid  analysis  of  conformal  antennas. 
This  hybrid  analysis  employs  an  integral  equation  formulation  for  dealing  with  antennas 
and  large  antenna  arrays  mounted  conformally  on  electrically  large,  perfectly-conducting 
convex  surfaces  with  a  uniform  material  coating  (which  can  act  as  a  radome  or  be  used  for 
impedance  matching,  or  both).  The  integral  equation  is  to  be  solved  numerically  via  the 
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moment  method  (MM);  however,  the  number  of  unknowns  to  be  solved  is  reduced  drastically 
by  a  special  choice  of  a  UTD  form  of  the  kernel  (Green’s  function)  in  the  integral  equation. 
This  analytical  UTD  Green’s  function  for  electric/magnetic  point  current  sources  on  an 
“arbitrary”  smooth  convex  conducting  surface  with  a  thin  uniform  material  coating  has 
been  completed  recently  and  constitutes  a  significant  result  which  is  to  be  used  in  this 
hybrid  approach.  The  general  UTD  solution  is  developed  from  generalization  of  the  canonical 
conducting  sphere  and  circular  cylinders  with  uniform  material  coating  which  was  reported  in 
the  last  year’s  JSEP  annual  report.  Presently,  a  numerical  computation  of  some  of  the  special 
UTD  transition  functions  is  being  conducted  so  that  this  analytical  UTD  Green’s  function 
can  be  computed  efficiently  for  use  in  the  hybrid  method.  It  is  noted  that  brute  force 
numerical  methods  may  become  intractable  for  predicting  the  behavior  of  large  scanning 
arrays  mounted  conformally  on  large  smooth  convex  conducting  bodies  with  or  without 
material  coating;  in  contrast  the  present  hybrid  method  is  expected  to  efficiently  analyze 
such  a  radiating  configuration  in  which  the  unknowns  to  be  solved  are  restricted  only  to  the 
fields  on  just  the  antennas  and  not  to  the  rest  of  the  large  convex  material  coated  platform. 
Demonstrations  of  the  utility  of  this  hybrid  approach  are  planned  for  the  next  phase  of  this 
study. 

Finally,  work  on  the  development  of  the  hybrid  procedure  for  analyzing  the  radia¬ 
tion/scattering  from  generic  but  realistic  complex  aircraft /missile  shapes  is  progressing  well. 
At  this  time,  the  required  UTD  solutions  which  will  comprise  the  special  Green’s  function 
for  the  fuselage  region  are  being  developed  in  the  unit  on  Diffraction  Studies.  These  UTD 
solutions  are  expected  to  be  completed  by  early  next  year.  They  will  subsequently  be  im¬ 
plemented  as  the  kernel  of  an  integral  equation  for  this  generic  radiating  configuration  to 
reduce  the  region  of  the  unknown  fields  to  be  solved  via  this  integral  equation  to  the  hori¬ 
zontal/vertical  stabilizers  and  fins,  since  the  main  fuselage  is  accounted  for  by  the  UTD.  The 
unknowns  over  the  regions  of  the  stabilizers  or  fins  will  then  be  solved  for  fairly  realistic  but 
generic  examples  of  aircraft /missile  using  the  moment  method  (MM).  The  MM  solution  will 
be  computed  via  an  iterative  approach  to  avoid  problems  associated  with  inversion  of  large 
MM  matrices.  In  the  future  phase  of  this  study,  a  super  hybrid  scheme  will  be  incorporated 
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to  include  the  effects  of  antennas,  antenna  windows,  gaps,  and  jet  engine  cavities.  Work  is 
in  progress  to  analyze  the  latter  components  by  techniques  best  suited  for  each  case. 

A  separate  yet  related  study  was  undertaken  to  substantially  speed  up  the  fast  multipole 
method  (FMM)  which  has  emerged  recently  as  a  very  useful  method  for  actually  speeding  up 
MM  computations  based  on  iterations.  This  acceleration  of  the  FMM  was  demonstrated  in 
the  2-D  case  by  using  a  hybrid  approach  in  the  sense  that  the  FMM  numerical  technique  was 
made  more  efficient  by  introducing  certain  windowed  operations  based  on  asymptotic  high 
frequency  concepts.  An  analysis  based  on  the  asymptotic  steepest  descent  method  provided 
a  rigorous  framework  for  incorporating  windows  into  the  FMM  scheme  thereby  substantially 
reducing  computation  times.  In  particular,  an  operation  count  0(N 3/2)  of  the  FMM  was 
reduced  to  0(N 4/3)  in  the  accelerated  FMM  procedure  developed  in  this  study.  Similar 
extensions  to  3-D  will  be  investigated  in  the  future  phases  of  this  study. 

2.  Design  and  Analysis  of  Active  Integrated  Antennas 

The  field  of  “wireless  communications”  has  experienced  a  tremendous  growth  in  the  last  few 
years  and  it  is  projected  that  this  growth  will  continue  in  the  foreseeable  future.  Some  of 
the  key  requirements  for  systems  in  this  field  such  as  cellular  phones  is  that  the  components 
need  to  be  low  cost,  small,  reliable,  etc.  Active  integrated  antennas  can  play  an  important 
role  because  they  satisfy  some  of  the  above  requirements.  The  integration  of  solid-state 
devices  (active)  and  antenna  functions  (passive)  onto  the  same  substrate  to  produce  compact 
radiating  structures  is  referred  to  as  an  active  integrated  antenna.  The  main  advantages  of 
such  structures  lie  in  their  low  cost,  simplicity  of  design  and  small  size.  Therefore,  they 
are  well-suited  for  mass  production  and  applications  in  radar,  satellite  and  communication 
systems  [l]-[4].  Some  of  the  goals  of  this  research  are  to  develop  analysis  and  design  methods 
that  are  accurate  and  efficient  as  well  as  designs  of  active  integrated  antennas  that  satisfy 
the  most  important  specifications  in  the  area  of  wireless  communications.  A  scheme  that  is 
currently  being  pursued  to  design/analyze  this  class  of  antennas  is  a  hybrid  combination  of  a 
full-wave  solver  (3-D  FDTD  code)  with  a  non-linear  circuit  simulator  (LIBRA,  developed  by 
EES  OF).  The  full  wave  solver  is  used  to  obtain  the  S-parameters  of  the  passive  components, 
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while  LIBRA  is  used  to  analyze  the  active  nonlinear  components.  The  full  wave  solver  will 
also  be  used  to  calculate  the  radiation  patterns  of  the  antennas. 

A  prototype  for  a  transistor-oscillator  active  antenna  was  developed  in  this  study  by 
integrating  a  three-terminal  device  (AT-42035,  medium  power  Silicon  Bipolar  Transistor, 
produced  by  HP-AVANTEK)  with  a  microstrip  patch  antenna  which  radiates  the  RF  power, 
and  acts  like  a  resonator  to  set  the  operating  frequency,  and  a  resonant  load  to  produce  high 
DC-to-RF  efficiencies.  Before  the  fabrication  process,  the  design  is  analyzed  and  verified 
with  an  improved  form  of  Kurokawa’s  generalized  analysis  [5]— [6] .  Note  that  Kurokawa’s 
analysis  is  usually  applied  without  taking  into  account  the  harmonics.  The  present  scheme 
being  considered  here  does  take  into  account  several  harmonics  and  therefore  yields  more 
accurate  designs. 

The  initial  design  can  be  done  by  the  “negative  resistance"  and/or  “reflection  coeffi¬ 
cient”  techniques  that  are  popular  in  the  design  of  oscillators.  We  define  ZA  and  T^  as  the 
impedance  and  the  reflection  coefficient  of  the  “active”  port,  and  similarly  Zl  and  Tj,  as  the 
impedance  and  the  reflection  coefficient  of  the  “load”  port  respectively.  For  the  reflection 
coefficient  technique,  the  oscillation  conditions  for  start-up  are: 

\Ta  II  rL  |>  1  (62) 

L(TA)  +  z(rL)  =  0.  (63) 

While  the  steady-state  conditions  can  be  written  as  follows: 

|Ta  II  Ft  |=  1  (64) 

Z(r^)  +  i(TL)  =  0  (65) 

For  the  negative  resistance  technique,  the  oscillation  conditions  for  start-up  are: 

Re(ZA)  +  Re{ZL)  <  0  (66) 

Im(ZA)  +  Im(ZL)  =  0  (67) 

and  the  steady-state  conditions  are  given  by: 

Re(ZA)  +  Re(ZL)  =  0  (68) 
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Im(ZA)  +  Im(ZL )  =  0  (69) 

Comparison  of  these  two  techniques  indicates  that  the  steady-state  conditions  for  both 
techniques  are  equivalent,  however  the  start-up  conditions  are  not,  and  the  negative  resis¬ 
tance  technique  is  much  better  for  design  purposes  [7].  The  reflection  coefficient  technique 
is  useful  for  analysis,  but  from  the  view  point  of  design,  the  oscillators  designed  based  on 
this  method  either  do  not  start  at  all  or  start  far  away  from  the  predicted  frequency.  On  the 
other  hand,  Equation  (66)  will  start  the  oscillation  and  then  the  non-linearities  will  bring  it 
into  steady-state.  Thus,  following  the  negative  resistance  technique,  the  following  steps  are 
important: 

•  An  NPN-active  biasing  circuit  is  precisely  modelled  in  Libra  and  a  non-linear  model  of 
AT-42035  is  biased  such  that  the  transistor  works  in  the  active  region.  Some  transistor 
parameters  are  compared  with  the  experimental  results  and  both  the  bias  network  and 
non-linear  model  of  the  transistor  are  tested. 

•  The  transistor  is  forced  to  become  potentially  unstable  at  the  frequency  of  oscillation 
by  adding  some  stubs  to  the  emitters  and  base. 

•  The  length  of  the  stubs  are  optimized  in  such  a  way  that  we  will  have  a  negative 
impedance  (“active”  port)  at  the  desired  frequency  and  the  emitters  have  a  balanced 
structure. 

The  output  port  of  the  active  device  is  directly  connected  to  a  microstrip  patch  antenna 
via  a  microstrip  line  (see  Figure  35)  and  the  design  of  the  load  should  satisfy  (66)  and  (67) 
simult  aneou  sly. 

The  input  impedance  of  the  antenna  at  resonance  is  chosen  to  be  exactly  the  same  as  the 
characteristic  impedance  of  the  microstrip  line  in  order  not  to  use  any  matching  networks. 
Therefore,  a  notch  is  introduced  and  the  width  of  the  notch  is  optimized  with  a  modified 
transmission  line  model.  The  design  is  verified,  and  both  the  time  and  frequency  domain 
characteristics  of  the  antenna  are  obtained  by  using  a  3-D  FDTD  code. 

As  mentioned  before,  we  are  using  an  improved  form  of  Kurokawa’s  generalized  analysis 
of  negative-resistance  oscillators  where  the  effect  of  harmonics,  which  has  not  been  included 
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Figure  36:  Schematic  of  an  improved  Kurokawa  analysis  which  includes  harmonics. 

before,  is  taken  into  consideration.  Besides  the  fundamental  signal,  the  harmonics  see  the 
main  load  (antenna  unit)  and  travel  back  and  forth  between  the  antenna  and  the  oscillator 
units  as  depicted  in  Figure  36.  Consequently,  this  is  a  very  realistic  representation  of  the 
actual  circuit.  In  Kurokawa’s  generalized  analysis,  a  stable  operation  point  for  a  current 
waveform  of  the  form  i(t)  =  A(t)  cos(u;t  +  4>(t))  is  given  by  ; 

BRA(A<,)dXL(u0)  aXA{A0)dRL(w0),  .  „ 

M~BA  gZ - BA - a^1  > 0  (7°) 

where  Za  =  Ra  AjXa  and  Ra  <  0  is  the  impedance  of  the  active  device,  and  ZjJ  =  Ri-\-jXi 

is  the  load  impedance.  In  our  improved  analysis  the  load  impedance  locus  will  not  change, 

however,  as  seen  from  Figure  37,  the  harmonics  (which  see  the  true  load  and  reflect  back) 

will  change  the  device  impedance  locus.  With  this  analysis  we  can  also  predict  the  current, 

voltage  and  power  levels  at  the  fundamental  frequency  and  harmonics. 

Two  different  active  antennas,  one  having  microstrip  line  open  circuits  (MLOC)  and 

the  other  having  microstrip  line  short  circuits  (MLSC)  at  the  emitters,  were  designed  (for 

operation  at  2.01GHz)  and  fabricated.  Both  designs  have  MLOC  at  the  base.  Figures  38 

and  39  show  the  measured  spectrum  of  the  signals  for  both  designs  at  the  desired  frequency. 

For  these  2  devices  the  center  frequency  is  seen  to  be  2.019  GHz  and  2.047  GHz,  respectively. 

The  final  step  was  to  measure  the  radiation  patterns  in  the  compact  range  at  The  Ohio  State 

University.  Figs  40  and  41  show  both  the  co-polarizations  and  the  cross-polarizations  in  the 
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Kurokawa's  Stability  Test 


Figure  37:  Locus  of  load  impedance  (antenna)  (o-o-o)  and  negative  of  active  device 
impedance  (x-x-x)  as  a  function  of  frequency.  The  load  impedance  is  plotted  for  var¬ 
ious  values  of  frequency  ranging  from  2  to  2.06  GHz.  The  negative  of  the  active  device 
impedance  is  plotted  for  various  values  of  collector  current  (AC)  ranging  from  0.925  to  27.95 
mA. 
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Figure  39:  Spectrum  of  measured  output  signal  for  MLSC  design 
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E  and  H-planes.  The  patterns  are  very  good  and  show  very  low  cross-polarized  components 
in  the  broadside  region  (—90°  to  90°). 

3.  Characteristic-Based  Finite  Volume  Time  Domain  (FVTD) 
Solver 


Since  the  late  80’s  there  has  been  a  continuing  effort  in  the  Computational  Fluid  Dynamics 
(CFD)  community  to  apply  their  numerical  algorithms  to  solve  Maxwell’s  equations  in  the 
time  domain.  Among  the  several  methods  developed  so  far,  a  Finite  Volume  Time  Domain 
(FVTD)  scheme,  originated  by  Shankar’s  group  [8,  9]  has  emerged  as  a  viable  scheme  for  its 
robustness  and  versatility.  In  the  Computational  Electromagnetics  Community,  the  Finite 
Difference  Time  Domain  (FDTD)  [10]  method,  originally  developed  by  Yee  [11],  has  on  the 
other  hand  already  attained  a  dominant  position  as  a  time  domain  Maxwell  solver  despite 
some  undesirable  numerical  characteristics.  Thus,  we  began  to  investigate  the  FVTD  method 
developed  in  [8,  9]  to  see  what  are  the  differences,  advantages  and  disadvantages  of  this 
method  over  the  well  known  FDTD.  There  are  very  few  papers  in  the  literature  where 
there  is  a  thorough  analysis  comparing  the  performance  of  these  two  methods  in  solving 
Maxwell’s  equations.  It  is  important  to  mention  that  there  are  other  FVTD  schemes  also 
being  developed  in  the  CFD  community  [12],  however,  at  the  present  time  we  are  only 
considering  one  of  them. 

The  FVTD  scheme  being  studied  here  is  based  on  numerical  algorithms  used  in  Compu¬ 
tational  Fluid  Dynamics.  Maxwell’s  equations  need  to  be  written  in  conservation  law  form 
which  is  commonly  done  in  the  CFD  community.  Then  the  quantities  we  are  dealing  with 
are  electric  fluxes  and  magnetic  fluxes  following  the  concepts  of  fluid  dynamics.  The  integral 
form  of  Maxwell’s  equations  can  be  written  in  the  form  of  a  conservation  law  as  follows 

iJIlQiV+JIs^-ii=IILSdV  <n> 


where  Q  is  a  vector  quantity  to  be  conserved,  J-  is  a  dyadic  flux  going  out  of  the  volume  V, 
and  S  is  the  source  vector.  For  Maxwell’s  equations, 
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where  X  is  the  identity  dyadic.  In  rectangular  coordinates,  X  —  xx  +  yy  +  zz.  Finally,  the 
source  vector  S  is  given  by 


The  differential  form  of  (71)  is  given  by 

jtQ-V  ■(?)  =  $  (74) 

Recall  that  for  linear,  isotropic,  inhomogeneous  media,  B  =  and  D  =  e(r)E. 

Equation  (74)  can  be  easily  modified  to  incorporate  general  curvilinear  coordinate  system 
which  enables  us  to  handle  complex  curvilinear  geometries  in  physical  space.  To  utilize  this 
body  conformal  feature  we  need  to  have  a  grid  generator  mapping  complex  geometries  in 
physical  space  to  a  simple  one  in  computational  space.  We  developed  and  tested  a  prelim¬ 
inary  version  of  a  computer  code  for  this  purpose.  Note  that  (71)  can  be  used  to  develop 
a  “unstructured”  grid  formulation,  which  has  some  advantages  over  the  “structured”  grid 
formulation  that  uses  a  curvilinear  coordinate  system. 

To  solve  Partial  Differential  Equations  (PDE)  we  need  a  time  integration  scheme.  Unlike 
the  FDTD,  which  uses  the  Leap-frog  method  in  time  integration,  the  FVTD  we  are  consid¬ 
ering  uses  a  modified  version  of  the  upwind  Lax-Wendroff  scheme  which  is  widely  used  in 
CFD  problems  [13] — [14] .  Here  the  scheme  is  based  on  characteristic  theory  of  signal  prop¬ 
agation.  For  hyperbolic  equations  like  Maxwell’s  equations,  upwind-based  schemes  can  be 
constructed  to  provide  the  right  amount  of  numerical  dissipation  to  achieve  both  stability 
and  accuracy  of  computation. 

Originally  these  upwind-based  methods  in  CFD  have  been  devised  to  handle  discontinu¬ 
ities  in  solutions  like  shock  phenomena.  In  Computational  Electromagnetics  these  features 
may  be  useful  in  the  study  of  wave  propagation  over  multilayer  materials  with  large  change  in 
the  dielectric  constants  from  one  layer  to  the  next.  The  conventional  FDTD  method  might 
have  difficulties  handling  this  situation.  A  second  advantage  of  FVTD  over  FDTD  is  its 
body  conformal  capability.  FDTD  uses  spatially  staggered  grids  that  can  only  approximate 
curvilinear  boundaries  in  stair  case  manner.  Though  there  have  been  some  efforts  in  the 
Computational  Electromagnetics  community  to  implement  the  FDTD  using  nonorthogonal 
and  unstructured  grids  [10],  their  achievements  at  the  present  time  cannot  be  compared 


to  those  used  in  FVTD  solvers  in  their  flexibility  and  robustness.  This  is  because  the  CFD 
community  has  a  long  tradition  of  treating  general  curvilinear  geometries  which  arise  in  fluid 
dynamics  [15]-  [18]  Of  course,  the  use  of  curvilinear  coordinates  adds  to  the  complexity  of 
the  solver.  A  third  possible  advantage  that  one  can  expect  from  FVTD  is  its  feasibility  in 
implementation  of  absorbing  boundary  conditions.  The  truncation  of  computational  space  is 
a  basic  dilemma  to  be  addressed  in  Computational  Electromagnetics.  Much  of  the  research 
effort  in  FDTD  has  been  devoted  to  the  development  of  appropriate  numerical  absorbing 
boundary  conditions  that  can  guarantee  minimum  reflections  in  the  outer  boundaries  of 
computational  spaces  at  a  reasonable  computational  cost  [19].  Since  FVTD  is  based  on 
characteristic  theory  of  signal  propagation,  it  has  the  direct  control  of  outgoing  waves.  With 
the  help  of  flexible  curvilinear  grid  systems  and  simple  absorbing  boundary  conditions,  one 
may  be  able  to  reduce  the  reflection  from  outer  boundaries  with  little  extra  computation. 
It  should  also  be  mentioned  that  the  FVTD  being  considered  here  is  suitable  for  wide  band 
applications. 

Some  of  the  reasons  that  have  made  the  EM  community  skeptical  about  this  method  are 
the  following.  First  of  all,  the  FVTD  is  fairly  sophisticated  both  in  theory  and  implemen¬ 
tation  and  inevitably  requires  large  computational  resources.  This  comes  from  the  fact  that 
original  algorithms  were  developed  to  treat  basically  nonlinear  discontinuous  situations  in 
fluid  dynamics.  But  today’s  demanding  electromagnetic  problems  can  include  nonlinear  and 
wide  band  applications  which  cannot  be  adequately  handled  by  existing  techniques.  One 
possible  application  of  this  method  is  in  Antenna  problems  especially  ones  with  complex 
curvilinear  geometries  and  high  contrast  dielectric  loading.  Currently  we  finished  a  simple 
one  dimensional  Maxwell  solver  in  homogeneous  media  and  are  developing  and  testing  a  2 
dimensional  solver  in  a  rectangular  grid  system.  As  an  example  of  the  1-D  solver,  Figure  42 
shows  the  propagation  of  Gaussian  pulses  in  homogeneous  space  with  different  Courant  num¬ 
bers  (CFL)  ranging  from  2  to  1.5.  Each  Gaussian  pulse  shows  no  significant  dispersion  after 
1000  time  steps.  We  are  currently  carrying  out  a  quantitative  comparison  of  FVTD  and 
FDTD  in  one  and  two  dimensions. 

In  developing  a  FVTD  code  for  antenna  applications  we  have  several  problems  to  address 
in  our  research.  The  first  one  is  the  extension  of  the  solver  to  2-D  and  3-D  in  general  curvi- 
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Figure  42:  Gaussian  pulse  propagation  for  different  Courant  numbers  (CFL)  where  CFL 
c  At/  A  x. 
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linear  coordinate  systems.  We  are  investigating  several  multi- dimension  algorithms  such  as 
the  Strang  Type  time  splitting  method.  Once  we  have  a  FVTD  solver  in  3-D,  we  will  need  a 
robust  3-D  grid  generator  which  may  take  much  more  effort.  Though  there  exists  commercial 
3-D  elliptic  generation  software  for  Aerospace  Engineering,  we  need  to  implement  our  own 
code  to  be  specifically  used  in  Antenna  Radiation  problems.  Another  important  point  is 
the  problem  of  excitation  which  is  an  integral  part  of  the  antenna  simulation  process.  Since 
the  nature  of  the  algorithms  in  FVTD  are  quite  different  from  the  conventional  FDTD,  we 
need  to  develop  methods  of  field  excitation  for  antenna  applications.  Note  that  most  of  the 
FVTD  applications  in  the  CFD  community  so  far  have  been  only  in  scattering  problems 
like  RCS  predictions  of  aircraft;  on  the  other  hand,  there  have  been  few  attempts  to  solve 
radiation  problems. 

Development  of  Elliptic  Grid  Generator  for  FVTD 

As  a  pre-processor  for  the  FVTD  method  discussed  here  there  is  a  need  for  body-fitting 
grid  generators  to  deal  with  curvilinear  geometries.  This  feature  is  very  important  for  the 
EM  solver  to  accurately  handle  curvilinear  boundaries  in  the  given  electromagnetic  problems. 
Note  that  there  have  been  extensive  research  on  various  aspects  of  numerical  grid  generation 
algorithms  for  the  last  two  decades.  Most  of  them  were  intended  for  use  in  Computational 
Fluid  Dynamics  where  curvilinear  boundaries  are  common  such  as  air-foil  problems.  Among 
them,  the  elliptic  grid  generation  algorithm  is  widely  used  in  many  engineering  problems 
for  its  versatility.  For  example,  for  the  two-dimensional  case,  in  the  elliptic  grid  scheme,  we 
need  to  find  the  transformations  £  =  £(x,y)  and  77  =  Tj(x,y)  to  satisfy  Poisson’s  equation, 
namely, 

V2  £  =  U  +  tw  =  P  (75) 

and 

V  V  =  Vxx  d"  1)yy  =  Q  (76) 

Inverting  the  above  equations,  yields, 

022®«  -  2012®£„  +  g22Xm  =  -022 ~  011  (77) 
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Figure  43:  Mapping  from  computational  to  physical  space. 

and 

922Vtt  —  2<7l23/£r)  +  922Vnr\  =  ~ ^22-P^  ~  1  Qj/rj  (78) 

where  gn  =  x\  +  y|,  £r12  =  x^xn  +  y ^  and  g22  =  X*  +  y*. 

Fig  43  shows  the  map  between  the  computational  and  physical  spaces. 

To  generate  elliptic  grids  we  need  to  solve  the  quasi-linear  elliptic  equations  (77)  and  (78) 
numerically.  Here  we  chose  the  Successive  Over  Relaxation(SOR)  method  to  solve  the  elliptic 
equations  because  it  is  relatively  easy  to  implement  it  while  the  computational  cost  is  still 
reasonable.  Figure  44  shows  the  computer  generated  grid  of  the  transformation  represented 
by  £(x,y)  and  i ](x,y). 

So  far  we  finished  the  implementation  of  a  preliminary  version  of  a  two  dimensional 
body-fitting  grid  generator  that  can  make  several  simple  canonical  geometries.  Fig  45  shows 
one  example  of  a  body-fitted  coordinate  system  in  2  dimensions.  Since  the  quality  and 
performance  of  a  body-fitted  grid  for  this  FVTD  method  can  be  estimated  properly  after 
we  incorporate  the  elliptic  body-fitted  grid  system  into  the  established  EM  solver  code,  we 
postponed  further  improvement  and  concentrated  our  efforts  in  the  implementation  of  the 
EM  solver.  For  a  complete  version  of  a  2-D  grid  generation  code  we  need  to  make  further 
improvements  such  as  refined  internal  grid  control  capability,  acceleration  of  computing  speed 
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Figure  44:  Grids  in  computational  and  physical  spaces. 

by  adaptive  parameter  decision,  user  friendly  interface  subroutines  for  complex  boundary 
definitions,  multizone  capability  and  extension  to  3-dimensional  grids. 

4.  Hybrid  Analysis  of  the  Scattering  by  Material  Bodies  and 
Antennas  with  Edges/ Junctions 

As  indicated  in  Section  III  of  the  Diffraction  Studies  Unit,  the  study  of  diffraction  by  material 
bodies  with  edges  and/or  junctions  is  important  due  to  the  widespread  use  of  composite 
materials  in  the  construction  of  airplanes,  satellites,  antenna  systems,  etc.  In  the  past  few 
years,  we  have  developed/improved  analytical  methods  to  study  this  class  of  diffraction 
phenomena  in  the  Diffraction  Studies  Unit.  Due  to  the  complex  nature  of  the  field  behavior 
near  material  surfaces  with  edges/junctions,  the  number  of  canonical  shapes  that  can  be 
studied  analytically  is  very  limited.  In  order  to  study  a  wider  class  of  canonical  shapes 
with  arbitrary  material  properties,  it  is  necessary  to  use  some  combination  of  analytical  and 
numerical  schemes.  The  scheme  chosen  here  was  the  Method  of  Moments/ Green’s  function 
approach  [20],  where  the  Green’s  function  of  a  canonical  shape  is  obtained  analytically  and 
the  rest  of  the  scatterer  is  handled  by  means  of  the  Method  of  Moments  (MM).  The  present 
study  was  done  for  a  material  cylinder  in  the  presence  of  a  wedge  shaped  object  as  depicted 
in  Figure  46.  An  special  case  of  the  material  cylinder  is  an  inhomogeneous  resistive  card 
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ure  45:  Curvilinear  coordinates  for  two-dimensional  body. 
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IMPEDANCE  WEDGE 


Figure  46:  Configuration  for  a  material  cylinder  in  the  presence  of  an  impedance  wedge. 

attached  to  the  edge  of  the  wedge  as  shown  in  Figure  47.  The  faces  of  the  wedge  satisfy 
impedance  (Leontovich)  boundary  conditions.  This  set  of  boundary  conditions  can  be  used 
to  model  coated  conductors,  rough  surfaces  or  surfaces  with  corrugations.  If  necessary, 
the  present  Green’s  function  can  be  modified  to  included  Generalized  Impedance  Boundary 
Conditions  (GIBC)  [21] 

The  development  of  the  Green’s  function  for  the  wedge  with  impedance  faces  is  discussed 
in  the  Diffraction  Studies  unit  and  will  not  be  repeated  here.  Using  this  special  Green’s 
function,  integral  equations  were  developed  for  the  geometries  depicted  in  Figures  46  and  47. 
This  integral  equations  are  then  solved  by  means  of  the  Method  of  Moments.  The  number 
of  unknowns  is  substantially  reduced  because  the  Green’s  function  automatically  takes  into 
account  the  presence  of  the  wedge.  However,  special  care  must  be  taken  in  solving  the 
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Figure  47:  Applications  for  an  inhomogeneous  resistive  card  to  reduce  the  scattering  from: 
(a)  a  PEC  semi-infinite  half-plane,  (b)  an  impedance  wedge. 


no 


integral  equation  for  Figure  47  for  the  T Ez  polarization  due  to  the  highly  singular  nature  of 
the  derivative  of  the  Green’s  function. 

An  important  application  of  this  study  is  the  analysis  of  the  radiation  characteristics  of 
antennas  mounted  on  a  finite  ground  plane.  The  radiated  field  will  in  general  illuminate 
the  edges  of  the  ground  plane  from  which  it  is  diffracted  over  a  wide  angular  sector.  This 
diffracted  field  can  then  serve  to  corrupt  the  radiation  pattern  of  the  antenna.  As  an  example, 
consider  a  microstrip  antenna  which  is  used  in  many  applications.  The  two-dimensional 
models  of  this  antenna  which  are  valid  in  the  E-  and  H-planes  are  shown  in  Figure  10  in 
the  Diffraction  Studies  unit.  One  way  the  edge  diffraction  can  be  decreased  is  to  attach  an 
inhomogeneous  resistive  card  to  the  edges  of  the  ground  plane.  The  present  analysis  can  be 
used  to  design  an  optimum  resistive  taper  because  it  is  a  very  efficient  scheme.  Note  that 
the  design  of  the  taper  is  a  synthesis  problem  which  in  general  has  to  be  solved  numerically 
due  to  the  nonlinear  nature  of  the  cost  function.  The  numerical  scheme  used  in  this  study  to 
optimize  the  resistive  card  is  the  Conjugate  Gradient  method  (CGM)  [22].  The  profile  of  the 
resistive  card  is  assumed  to  be  exponential,  quadratic  or  some  other  function.  This  reduces 
the  number  of  unknowns  to  two  or  three  variables.  However,  if  no  function  is  assumed  for 
the  taper,  the  number  of  unknowns  is  larger.  We  are  currently  planning  to  use  the  Genetic 
Algorithm  [23]  which  is  more  efficient  than  the  CGM  when  there  are  a  large  number  of 
variables  that  are  being  optimized. 

Although  the  MM/Green’s  function  method  developed  here  is  applicable  to  a  structure 
with  a  single  edge,  it  can  be  easily  adapted  to  study  structures  with  two  or  more  edges. 
This  is  accomplished  by  obtaining  a  numerical  diffraction  coefficient  for  a  single  edge  and 
then  using  this  diffraction  coefficient  to  diffract  from  various  edges.  Consider  the  single¬ 
edge  configuration  shown  in  Figure  48(a);  where  the  total  radiated  field  consists  of  the  two 
field  components  Ugo,  which  is  the  geometrical  optics  field  that  would  be  radiated  by  an 
antenna  that  is  mounted  on  an  infinite  substrate  (shown  in  Figure  48(b)),  and  U^,  which 
is  the  field  that  diffracts  from  the  first  edge.  The  pattern  for  an  antenna  mounted  on  a 
structure  with  a  single  edge,  and  whatever  material  loading  there  may  be,  is  generated  using 
the  MM/Green’s  function  technique.  Then,  the  pattern  for  the  same  antenna,  but  this  time 
mounted  on  an  infinite  substrate,  is  then  generated  by  using  the  wedge  Green’s  function.  By 
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Figure  48:  Development  of  the  hybrid  MM/GTD  technique.  The  single-edge  geometry 
consists  of  the  fields  Ugo  and  (a).  The  field  Ugo  is  determined  from  the  infinite  geometry 
(b).  Combining  the  field  term  Ud2  to  the  single-edge  geometry  in  (a)  then  gives  the  pattern 
for  the  double-edge  geometry  (c). 
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subtracting  the  pattern  of  the  infinite  structure  from  the  pattern  for  the  single-edge  geometry, 
the  diffracted  field  term  is  obtained.  By  shifting  the  phase  and  angle  reference,  this  term 
can  then  be  used  to  obtain  the  field  that  diffracts  from  the  second  edge  Ud3,  which  is  then 
superimposed  onto  the  pattern  for  the  single-edge  geometry  (Figure  48(a))  to  give  the  pattern 
for  the  double-edge  geometry  (Figure  48(c)).  This  method  will  be  referred  to  as  the  hybrid 
MM/GTD  method.  Although  the  pattern  which  has  been  generated  by  this  method  does 
not  take  into  account  multiple  diffraction  that  can  occur  between  the  two  edges,  it  will  be 
shown  that  the  single  diffraction  model  is  a  very  accurate  one  and  can  be  extremely  useful  for 
developing  methods  to  control  the  edge  diffraction.  The  configuration  for  an  actual  antenna 
that  has  been  examined  using  the  present  approach  is  depicted  in  Figure  49.  The  original 
antenna  had  an  all  ground  plane  and  it  was  the  goal  of  this  design  to  replace  a  section  of  that 
ground  plane  by  tapered  resistive  cards.  The  two  dimensional  model  that  is  used  for  this 
finite  configuration  is  shown  in  Figure  50.  This  2-D  model  is  valid  in  the  E-  and  H-planes 
where  the  patch  is  replaced  by  the  appropriate  z-directed  current  source.  Note  that  the 
grounded  dielectric  substrate  is  modeled  by  an  equivalent  impedance  surface.  One  of  the 
advantages  of  using  the  wedge  Green’s  function  in  this  analysis  is  that  the  different  field 
components  can  be  explicitly  separated  and  studied  individually.  This  can  be  a  valuable 
tool  for  providing  physical  insight  into  the  various  mechanisms  involved  in  the  radiation  and 
scattering  from  structures.  In  determining  the  radiation  pattern  of  the  patch  in  the  presence 
of  the  edge,  there  are  three  field  components  that  should  be  considered,  which  are  shown 
in  Figure  51.  Usually,  the  direct  space  wave  is  the  desired  field  and  the  other  two,  namely, 
diffracted  space  and  surfaces  waves  are  not  desired  because  they  distort  the  direct  space 
wave.  Radiation  patterns  will  be  shown  for  the  antenna  configuration  depicted  in  Figure  49, 
where  the  entire  ground  plane  is  a  PEC  and  when  a  portion  of  the  ground  plane  is  replaced 
by  an  optimized  resistive  card  with  a  quadratic  taper.  For  comparison,  measured  results 
as  well  as  MM  reference  solution  are  also  included.  This  MM  reference  solution  utilizes  a 
free-space  Green’s  function. 

A  plot  of  the  radiation  pattern  for  the  all  PEC  ground  plane  configuration  comparing 
the  hybrid  MM/GTD  and  MM  solutions  for  the  2-D  model  with  measured  data  is  shown  in 
Figure  52.  The  agreement  between  the  hybrid  MM/GTD  and  MM  solutions  is  very  good, 
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Figure  49:  Three-dimensional  configuration  for  a  rectangular  microstrip  patch  radiator 
mounted  on  a  grounded,  dielectric  slab. 


(a) 

- x 


(b) 


Figure  50:  2-D  configuration  for  the  finite  antenna  of  Figure  49.  Dielectric  substrate  can  be 
mounted  on  an  entirely  PEC  ground  plane  (a),  or  on  a  ground  plane  where  part  of  the  PEC 
is  replaced  by  an  R-card  (b). 
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Figure  51:  Various  field  components  for  a  radiating  patch  on  a  substrate. 

and  the  slight  discrepancies  can  be  mostly  explained  by  the  lack  of  multiple  diffraction  terms 
in  the  hybrid  method.  The  patterns  generated  by  the  2-D  model  compare  reasonably  well 
with  the  measured  data,  predicting  the  same  characteristic  shape  and  ripple  in  the  mainbeam 
as  well  as  the  back  region,  with  some  noticeable  level  difference  there.  This  is  to  be  expected 
since  there  are  some  difficulties  in  measuring  the  radiation  pattern  in  the  back  region,  mainly 
due  to  the  presence  of  the  mounting  structure  and  feed  network  for  the  antenna.  There  is  also 
some  asymmetry  in  the  measured  pattern,  which  indicates  there  could  be  some  measurement 
error  as  well.  A  plot  showing  the  same  comparisons  is  next  made  for  the  case  where  the 
R-card  with  the  quadratic  profile  replaces  a  section  of  the  ground  plane  (see  Figure  53).  In 
this  case,  the  agreement  between  the  hybrid  and  MM  methods  is  noticeably  better.  This  is 
most  likely  due  to  the  reduction  in  the  amount  of  multiple  diffraction  because  of  the  presence 
of  the  R-cards.  The  R-cards  reduce  the  end  point  diffraction  with  the  obvious  benefit  seen 
by  the  reduction  of  the  ripple  in  the  mainbeam,  which  is  shown  in  both  the  2-D  model  and 
the  measurements,  and  it  also  decreases  the  amount  of  scattering  in  the  back  region.  The 
improvement  of  the  antenna  performance  in  the  mainbeam  region  is  more  clearly  shown  by 
the  phase  plot  of  Figure  54  for  the  measured  pattern.  The  diffraction  from  the  edges  of 
the  PEC  ground  plane  introduces  an  oscillating  phase  variation  that  indicates  an  unstable 
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Figure  52:  TE2  radiation  patterns  for  the  finite  configuration  of  Figure  50(a)  (with  L  =  .3A. 
L\  —  .67A,  L 3  =  2.54A,  and  t  =  .0184A),  comparing  the  hybrid  MM/GTD  and  MM  solutions 


for  the  2-D  model  with  measured  data  for  an  all  PEC  ground  plane. 
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Figure  53:  TEZ  radiation  patterns  for  the  finite  configuration  of  Figure  50(b)  (with  L  =  ,3A, 
Li  =  .67A,  L2  =  1A,  L3  =  2.54A,  and  t  =  .0184A),  comparing  hybrid  MM/GTD  and  MM 
solutions  for  the  2-D  model  with  measured  data  for  an  R-card  with  a  quadratic  profile. 


Figure  54:  TEZ  measured  phase  responses  for  the  configurations  of  Figures  50(a)  and  (b) 
(with  L  =  .3A,  fa  =  .67A,  £3  =  2.54A,  and  t  —  .0184A),  comparing  the  effect  of  replacing 
the  PEC  with  an  R-card. 

phase  center  for  the  antenna,  whereas  the  R-card  provides  for  a  very  flat  phase  characteristic 
which  indicates  a  more  stable  phase  center.  Similar  patterns  can  also  be  calculated  for  the 
H-plane,  however,  they  will  not  be  shown  here. 

5.  Conformal  Antennas 

A  hybrid  MM-UTD  approach  was  developed  previously  under  JSEP  to  efficiently  analyze 
large  microstrip  or  slot  arrays  on  planar  conducting  surfaces  with  a  plane  material  coating. 
This  was  made  possible  by  including  an  asymptotic  or  UTD  form  of  the  special  grounded 
material  slab  Green’s  function  which  occurs  as  the  kernel  of  the  governing  integral  equation 
for  solving  the  unknown  currents  on  the  microstrip  patches  or  the  fields  in  the  slots.  The 
unknowns  are  then  solved  by  the  MM  technique.  Since  the  UTD  form  of  the  Green’s  function 
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is  in  closed  form  as  opposed  to  the  exact  conventional  Sommerfeld  integral  form  which  is 
poorly  convergent,  the  hybrid  MM-UTD  makes  the  analysis/synthesis  of  truly  electrically 
large  arrays  tractable,  as  opposed  to  the  use  of  a  conventional  brute  force  numerical  MM 
solution  with  the  exact  Sommerfeld  integral  form  of  the  special  Green’s  function  to  do  the 
same.  It  is  noted  that  infinite  periodic  array  theory  which  is  often  employed  to  study  such 
arrays  is  not  strictly  applicable  because  practical  arrays  are  of  finite  extent  and  the  effects 
of  finiteness  of  the  array  are  included  only  crudely  in  such  a  infinite  array  model.  On  the 
other  hand,  the  present  hybrid  approach  rigorously  accounts  for  the  finiteness  of  the  array 
and  predicts  the  array  edge  effects  accurately.  Work  is  now  in  progress  to  employ  this  UTD 
special  Green’s  function,  whose  usefulness  has  already  been  demonstrated  for  analysis  of 
one-dimensional  arrays,  to  next  deal  with  two-dimensional  arrays  containing  feed  networks, 
etc.,  to  demonstrate  the  versatility,  efficiency  and  usefulness  of  this  hybrid  approach  for 
analysis/design  of  large  finite  arrays. 

Yet  another  important  hybrid  MM-UTD  solution  is  being  developed  to  deal  with  arrays 
on  non-planar  convex  surfaces  of  variable  curvature.  In  particular,  as  a  first  step,  a  UTD 
form  of  the  special  Green’s  function  has  been  developed  for  a  fully  3-D  material  coated 
perfectly  conducting  convex  surface,  which  for  the  first  time  would  predict  the  performance 
of  antennas  in  or  on  such  configurations  and  show  explicitly  in  relatively  closed  form,  the 
effects  on  fields  of  such  antennas  arising  from  surface  curvature  and  the  presence  of  material 
coating.  It  is  noted  that  infinite  periodic  array  theory  is  not  even  applicable  to  such  curved 
surface  configurations  because  arbitrary  curvature  destroys  periodicity  in  the  fields;  further¬ 
more,  it  is  note  possible  to  obtain  an  exact  analytical  form  of  the  special  Green’s  function 
for  this  case.  Hence,  the  importance  of  developing  a  special  UTD  based  special  Green’s 
function  for  an  arbitrary  smooth  convex  surface  with  uniform  material  coating  can  hardly 
be  overemphasized. 

In  the  recent  months,  the  UTD  transition  functions  involving  modified  Fock  type  integrals 
that  appear  in  the  UTD  form  of  the  special  Green’s  function,  which  has  been  developed  for 
electric  or  magnetic  point  sources  on  an  arbitrary  perfectly  conducting  convex  surface  with  a 
uniform  material  coating,  are  being  computed  numerically  in  an  efficient  manner.  While  some 
past  researchers  have  evaluated  some  of  the  these  types  of  Fock  integrals,  their  procedure  is 


120 


extremely  time  consuming.  Several  useful  procedures  have  been  incorporated  to  significantly 
speed  up  the  evaluation  of  these  integrals,  and  to  also  enhance  the  accuracy  of  computations. 
It  is  imperative  that  these  special  UTD  transition  functions  be  computed  rapidly  to  maintain 
the  overall  efficiency  and  utility  that  are  characteristic  of  the  UTD  method  in  general.  There 
are  four  basic  types  of  UTD  Fock  integral  transition  functions  which  are  required  in  the 
analysis  of  radiation  and  mutual  coupling  associated  with  antennas  on  perfectly-conducting 
convex  surface  with  uniform  material  coating.  Unlike  the  Fock  integral  based  UTD  transition 
functions  which  are  universal  functions  of  only  a  single  parameter  for  the  perfectly  conducting 
convex  surface  (and  which  hence  can  be  tabulated),  the  ones  for  the  material  coated  convex 
surface  depend  on  at  least  two  parameters.  Furthermore,  their  behavior  is  more  complicated 
than  for  the  uncoated  case.  Three  out  of  these  four  basic  Fock  type  UTD  transition  functions 
have  been  coded  efficiently  using  some  useful  procedures  to  sum  the  “tails”  of  these  infinite 
Fock  integrals  in  closed  form.  The  remaining  fourth  basic  Fock  type  function  for  the  coated 
case  is  currently  under  investigation  to  study  its  properties  carefully  so  that  an  intelligent 
choice  of  procedures  can  then  be  developed  to  evaluate  this  Fock  type  integral  efficiently  as 
well.  It  is  expected  that  the  efficient  numerical  evaluation  of  these  special  UTD  transition 
functions  will  be  completed  by  the  end  of  this  year  or  so.  Subsequently,  it  is  proposed  to 
employ  this  UTD  solution  in  a  hybrid  MM-UTD  formulation  to  demonstrate  its  versatility 
and  usefulness  in  the  analysis/design  of  antenna  arrays  that  truly  conform  to  a  general 
smooth  doubly  curved  convex  surface  (e.g.,  the  fuselage  of  an  aircraft  or  missile,  or  the  top 
of  an  automobile,  etc.).  The  antenna  elements  can  be  relatively  arbitrary  (slots,  monopoles, 
microstrip  patches,  etc.)  in  this  hybrid  MM-UTD  approach.  It  is  important  to  note  that  no 
such  solution  presently  exists,  other  than  this  proposed  hybrid  MM-UTD  approach,  for  the 
analysis/design  of  truly  general  convexly  curved,  large  conformal  antenna  arrays. 

6.  High-Frequency  Asymptotic  Acceleration  of  the  Fast  Multi¬ 
pole  Method 

The  fast  multipole  method  (FMM)  is  a  fairly  recent  development  [24]  which  provides  a 
means  of  accelerating  iterative  algorithms  for  solving  large-scale  integral  equation  problems. 
In  our  research,  the  FMM  has  been  improved  so  that  it  is  substantially  faster  as  will  be 
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explained  below.  This  FMM  approach  may  be  applied  to  electromagnetic  integral  equation 
problems  where  the  unknown  currents  are  expanded  into  N  subsectional  basis  elements.  It¬ 
erative  algorithms,  such  as  the  conjugate  gradient  method,  solve  for  the  unknown  currents 
by  repeatedly  multiplying  the  N  x  N  system  matrix  times  the  0(N )  current  vector.  This 
matrix-vector  product  computation  has  an  operational  count  of  0(N2),  and  many  iterations 
may  be  required  for  convergence.  The  basic  FMM  accelerates  the  matrix-vector  multipli¬ 
cation  by  dividing  the  N  basis  elements  into  M  localized  groups,  and  by  expressing  the 
matrix  interactions  between  separated  groups  using  multipole  expansions.  The  matrix  inter¬ 
actions  between  basis  elements  inside  the  same  group  and  with  elements  of  adjacent  groups 
are  expressed  in  the  conventional  manner,  so  singularities  associated  with  certain  integral 
equation  formulations  may  be  handled  using  existing  methods.  By  optimally  choosing  M 
to  minimize  the  number  of  computer  operations,  the  FMM  reduces  the  operational  count  of 
the  matrix-vector  product  from  0(N2)  to  0(N 3/2). 

The  matrix  interactions  between  separated  groups  in  the  FMM  are  expressed  via  plane 
wave  translations  from  the  source  group  to  the  receiving  group.  In  the  basic  FMM,  plane 
waves  propagate  in  all  directions  and  all  plane  waves  are  translated  from  the  source  to 
the  receiver.  However,  from  a  high-frequency  asymptotic  point  of  view,  only  plane  waves 
propagating  close  to  the  direction  from  source  to  receiver  should  contribute  strongly.  A 
windowed  version  of  the  FMM  has  been  developed  which  filters  out  plane  waves  propagating 
at  wide  angles  with  respect  to  the  line  connecting  two  separated  groups.  The  windowing 
filter  is  derived  from  an  asymptotic  evaluation  of  the  plane  wave  translation  operator.  The 
resulting  reduction  in  plane  waves  further  accelerates  the  FMM  to  an  0(N 4/3)  operational 
count  without  increasing  the  complexity  of  the  implementation.  Figure  55  shows  the  CPU 
times  for  a  single  matrix-vector  multiplication  found  using  brute  force,  the  basic  FMM, 
and  the  new  asymptotically  windowed  FMM.  The  solid  lines  are  actual  times  on  a  Sun 
SparcStation  IPX,  and  the  dashed  lines  are  estimated  projections. 


122 


Figure  55:  CPU  time  vs.  number  of  unknowns  N  for  a  single  matrix-vector  multiplication. 
Solid  lines  are  computed  data  points  and  dashed  lines  are  estimated  projections.  Times  are 
for  a  Sun  SPARCstation  IPX  workstation. 
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7.  Hybrid  High-Frequency  Asymptotic/Numerical  Methods  for 
Electrically  Large  Radiating  Objects 

An  important  EM  phenomenon  that  needs  to  be  predicted  efficiently  and  accurately  is  the  ra¬ 
diation/scattering  of  EM  waves  by  electrically  large  realistic  aircraft/missile  geometry  which 
obviously  falls  in  the  category  of  complex  shapes.  Such  a  complex  realistic  aircraft  /missile 
configuration  typically  consists  of  a  large  smooth  and  generally  convex  fuselage  except  possi¬ 
bly  for  an  end  cap  in  the  nose  region  to  use  a  forward  looking  antenna  inside  a  radome,  and 
has  large  control  surfaces  (wings,  fins,  etc.)  as  well  as  other  small  and  large  perturbations 
such  as  antennas,  antenna  windows,  jet  inlet  cavities,  gaps  etc.  A  brute  force  numerical 
EM  solution  of  such  an  electrically  large,  complex  radiating/scattering  object  can  become 
tremendously  inefficient  and  even  intractable.  On  the  other  hand,  a  purely  high  frequency 
based  EM  solution  may  not  be  able  to  handle  electrically  small  but  nonetheless  potentially 
important  contributors  to  the  overall  radiation/scattering,  and  it  also  may  not  be  able  to 
handle  some  electrically  large  portions  for  which  the  pertinent  asymptotic  high  frequency 
solutions  (e.g.,  the  UTD)  are  not  available  at  the  present  time.  Clearly,  therefore,  a  hybrid 
MM-UTD  approach  appears  to  provide  an  attractive  scheme  for  analyzing  this  complex  EM 
configuration,  where  MM  refers  to  the  numerical  moment  method  (MM)  for  solving  the  gov¬ 
erning  integral  equation,  and  the  kernel  of  the  integral  equation  (pertinent  Green’s  function) 
is  to  be  approximated  in  this  hybrid  approach  by  UTD  to  a  priori  account  for  the  basic 
convex  aircraft /missile  fuselage  shape  leaving  the  unknown  induced  surface  currents  to  be 
found  only  over  the  remaining  part  of  the  aircraft /missile  geometry.  Such  a  hybrid  MM- 
UTD  approach  is  far  more  efficient  than  a  brute  force  numerical  MM  solution  of  the  integral 
equation  for  the  whole  geometry  or  a  numerical  solution  based  on  the  partial  differential 
equations  of  EM,  because  the  number  of  unknowns  to  be  solved  in  the  hybrid  approach  are 
drastically  reduced  in  comparison  to  those  that  need  to  be  found  in  the  brute  force  numerical 
solution  methods.  However,  one  notes  that  even  with  this  hybrid  procedure,  it  is  possible  to 
still  have  a  large  number  of  unknowns,  even  though  they  may  be  far  less  than  the  number 
of  unknowns  required  to  be  solved  by  conventional  numerical  based  methods.  It  is  for  this 
reason  one  must  not  only  find  the  best  choice  of  the  kernel  (Green’s  function),  which  is  to  be 
represented  asymptotically  via  a  UTD  approximation,  in  order  to  a  priori  limit  the  region 
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of  the  unknowns  to  be  solved  using  the  MM  technique,  but  one  must  in  addition  also  find 
ways  to  speed  up  the  MM  portion  of  the  computations.  In  particular,  the  latter  implies  an 
efficient  choice  of  basis  or  expansions  functions  for  the  unknown  currents  and  an  iterative 
solution  of  the  MM  equation  rather  than  the  conventional  MM  matrix  inversion  method  of 
solution.  Presently,  it  is  planned  that  the  UTD  solutions  being  developed  for  the  smooth 
convex  surface  (with  or  without  an  end  cap)  for  arbitrary  source  and  field  point  locations  as 
described  in  Part  IV  of  the  section  on  Diffraction  Studies,  will  be  incorporated  into  a  version 
of  our  existing  aircraft  geometry  code  which  is  also  simultaneously  undergoing  additional 
improvement  in  its  modeling  capabilities.  It  is  expected  that  the  relevant  UTD  results  will 
be  incorporated  by  mid- 1996  into  this  aircraft  geometry  code.  This  geometry  code  will  also 
provide  a  discretization  of  the  surface  geometry  as  required  for  the  MM  part  of  this  hybrid 
scheme.  Once  this  basic  hybrid  MM-UTD  scheme  is  demonstrated  for  a  more  sophisticated 
generic  aircraft  shape  than  one  consisting  of  just  plates  and  circular  cylinders,  it  will  later  be 
extended  through  a  super-hybrid  scheme  to  include  additional  effects  such  as  those  resulting 
from  the  presence  of  an  inlet  cavity  or  antennas,  etc.  It  is  noted  that  separate  procedures 
are  being  developed  to  deal  with  inlets,  gaps/cracks  and  antennas,  etc.  This  super  hybrid 
scheme  will  be  highly  useful  once  it  is  fully  developed,  for  it  is  expected  that  this  may  ac¬ 
tually  be  one  of  the  most  tractable  approaches  for  analyzing  the  radiation/scattering  from 
generic  and  also  actual  aircraft /missile  shapes. 
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